Fisher线性判别

在理解Fisher线性分类的参考文件的代码基础上（matlab代码），改用python代码完成Fisher判别的推导。重点理解“群内离散度”（样本类内离散矩阵）、“群间离散度”（总类内离散矩阵）的概念和几何意义。

代码实现

matlab代码

clc
clear
data=xlsread('Iris.csv');
Iris1=data(1:50,1:4);
Iris2=data(51:100,1:4);
Iris3=data(101:150,1:4);
%类均值向量
m1 = mean(Iris1);
m2 = mean(Iris2);
m3 = mean(Iris3);
 
%各类内离散度矩阵
s1 = zeros(4);
s2 = zeros(4);
s3 = zeros(4);
for i=1:1:30
    s1 = s1 + (Iris1(i,:) - m1)'*(Iris1(i,:) - m1);
end
for i=1:1:30
    s2 = s2 + (Iris2(i,:) - m2)'*(Iris2(i,:) - m2);
end
for i=1:1:30
    s3 = s3 + (Iris3(i,:) - m3)'*(Iris3(i,:) - m3);
end
 
%总类内离散矩阵
sw12 = s1 + s2;
sw13 = s1 + s3;
sw23 = s2 + s3;
%投影方向
w12 = ((sw12^-1)*(m1 - m2)')';
w13 = ((sw13^-1)*(m1 - m3)')';
w23 = ((sw23^-1)*(m2 - m3)')';
%判别函数以及阈值T（即w0）
T12 = -0.5 * (m1 + m2)*inv(sw12)*(m1 - m2)';
T13 = -0.5 * (m1 + m3)*inv(sw13)*(m1 - m3)';
T23 = -0.5 * (m2 + m3)*inv(sw23)*(m2 - m3)';

%% 请补充输出判别函数



%% 请给下面代码添加必要注释 
kind1 = 0;
kind2 = 0;
kind3 = 0;
newiris1=[];
newiris2=[];
newiris3=[];
for i=31:50
    x = Iris1(i,:);
    g12 = w12 * x' + T12;
    g13 = w13 * x' + T13;
    g23 = w23 * x' + T23;
    if((g12 > 0)&(g13 > 0))
       newiris1=[newiris1;x];
        kind1=kind1+1;
    elseif((g12 < 0)&(g23 > 0))
         newiris2=[newiris2;x];
    elseif((g13 < 0)&(g23 < 0))
          newiris3=[newiris3;x];
    end
end    
for i=31:50
    x = Iris2(i,:);
    g12 = w12 * x' + T12;
    g13 = w13 * x' + T13;
    g23 = w23 * x' + T23;
    if((g12 > 0)&(g13 > 0))
      newiris1=[newiris1;x];
    elseif((g12 < 0)&(g23 > 0))
            kind2=kind2+1;
           newiris2=[newiris2;x];
    elseif((g13 < 0)&(g23 < 0))
            newiris3=[newiris3;x];
    end
end
for i=31:50
    x = Iris3(i,:);
    g12 = w12 * x' + T12;
    g13 = w13 * x' + T13;
    g23 = w23 * x' + T23;
    if((g12 > 0)&(g13 > 0))
       newiris1=[newiris1;x];
    elseif((g12 < 0)&(g23 > 0))
           newiris2=[newiris2;x];
    elseif((g13 < 0)&(g23 < 0))
            kind3=kind3+1;
            newiris3=[newiris3;x];
    end
end

correct=(kind1+kind2+kind3)/60;
fprintf('\n综合正确率：%.2f%%\n\n',correct* 100);

python代码

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt 
import seaborn as sns

path=r'Iris.csv'
df = pd.read_csv(path, header=0)
Iris1=df.values[0:50,0:4]
Iris2=df.values[50:100,0:4]
Iris3=df.values[100:150,0:4]
m1=np.mean(Iris1,axis=0)
m2=np.mean(Iris2,axis=0)
m3=np.mean(Iris3,axis=0)
s1=np.zeros((4,4))
s2=np.zeros((4,4))
s3=np.zeros((4,4))
for i in range(0,30,1):
    a=Iris1[i,:]-m1
    a=np.array([a])
    b=a.T
    s1=s1+np.dot(b,a)    
for i in range(0,30,1):
    c=Iris2[i,:]-m2
    c=np.array([c])
    d=c.T
    s2=s2+np.dot(d,c) 
    #s2=s2+np.dot((Iris2[i,:]-m2).T,(Iris2[i,:]-m2))
for i in range(0,30,1):
    a=Iris3[i,:]-m3
    a=np.array([a])
    b=a.T
    s3=s3+np.dot(b,a) 
sw12=s1+s2
sw13=s1+s3
sw23=s2+s3
#投影方向
a=np.array([m1-m2])
sw12=np.array(sw12,dtype='float')
sw13=np.array(sw13,dtype='float')
sw23=np.array(sw23,dtype='float')
#判别函数以及T
#需要先将m1-m2转化成矩阵才能进行求其转置矩阵
a=m1-m2
a=np.array([a])
a=a.T
b=m1-m3
b=np.array([b])
b=b.T
c=m2-m3
c=np.array([c])
c=c.T
w12=(np.dot(np.linalg.inv(sw12),a)).T
w13=(np.dot(np.linalg.inv(sw13),b)).T
w23=(np.dot(np.linalg.inv(sw23),c)).T
#print(m1+m2) #1x4维度  invsw12 4x4维度  m1-m2 4x1维度
T12=-0.5*(np.dot(np.dot((m1+m2),np.linalg.inv(sw12)),a))
T13=-0.5*(np.dot(np.dot((m1+m3),np.linalg.inv(sw13)),b))
T23=-0.5*(np.dot(np.dot((m2+m3),np.linalg.inv(sw23)),c))
kind1=0
kind2=0
kind3=0
newiris1=[]
newiris2=[]
newiris3=[]
for i in range(30,49):
    x=Iris1[i,:]
    x=np.array([x])
    g12=np.dot(w12,x.T)+T12
    g13=np.dot(w13,x.T)+T13
    g23=np.dot(w23,x.T)+T23
    if g12>0 and g13>0:
        newiris1.extend(x)
        kind1=kind1+1
    elif g12<0 and g23>0:
        newiris2.extend(x)
    elif g13<0 and g23<0 :
        newiris3.extend(x)
#print(newiris1)
for i in range(30,49):
    x=Iris2[i,:]
    x=np.array([x])
    g12=np.dot(w12,x.T)+T12
    g13=np.dot(w13,x.T)+T13
    g23=np.dot(w23,x.T)+T23
    if g12>0 and g13>0:
        newiris1.extend(x)
    elif g12<0 and g23>0:
 
        newiris2.extend(x)
        kind2=kind2+1
    elif g13<0 and g23<0 :
        newiris3.extend(x)
for i in range(30,50):
    x=Iris3[i,:]
    x=np.array([x])
    g12=np.dot(w12,x.T)+T12
    g13=np.dot(w13,x.T)+T13
    g23=np.dot(w23,x.T)+T23
    if g12>0 and g13>0:
        newiris1.extend(x)
    elif g12<0 and g23>0:     
        newiris2.extend(x)
    elif g13<0 and g23<0 :
        newiris3.extend(x)
        kind3=kind3+1
correct=(kind1+kind2+kind3)/60
print("样本类内离散度矩阵S1：",s1,'\n')
print("样本类内离散度矩阵S2：",s2,'\n')
print("样本类内离散度矩阵S3：",s3,'\n')
print('------------------------------------------------------------------------------')
print("总体类内离散度矩阵Sw12：",sw12,'\n')
print("总体类内离散度矩阵Sw13：",sw13,'\n')
print("总体类内离散度矩阵Sw23：",sw23,'\n')
print('------------------------------------------------------------------------------')
print('判断出来的综合正确率：',correct*100,'%')

………

---

博客刚开，多多包涵，以后会尽量完善的
I’m Ray.I’m ok.