fisher.m

这是一各关于线性判别函数分类的MATLAB程序

clear
clc
close all;
%m1,m2均为10个样本
% m1=[-0.4,0.58,0.089;-0.31,0.27,-0.04;0.38,0.055,-0.035;-0.15,0.53,0.011;-0.35,0.47,0.034;0.17,0.69,0.1;-0.011,0.55,-0.18];
% m2=[0.83,1.6,-0.014;1.1,1.6,0.48;-0.44,-0.41,0.32;0.047,-0.45,1.4;0.28,0.35,3.1;-0.39,-0.48,0.11;0.34,-0.079,0.14];
% a) 对类别m1,m2,利用fisher原则计算投影方向w
m1=[0 0 0; 1 0 0; 1 0 2; 1 1 0; 0.5 0 0.5;1 0.5 0;1 0 0.5];
m2=[0 0 1; 0 1 0; 0 1 1; 1 2 1;0 0.5 1;0 1 0.5; 0.5 1 2];
u1=mean(m1);%求均值
u2=mean(m2);
%计算类内散度Si和总类内散度Sw
S1=(m1(1,:)-u1)'*(m1(1,:)-u1);
S2=(m2(1,:)-u2)'*(m2(1,:)-u2);
for i=2:7
    S1=S1+(m1(i,:)-u1)'*(m1(i,:)-u1);
    S2=S2+(m2(i,:)-u2)'*(m2(i,:)-u2);
end
Sw=S1+S2;
%变换向量
w=inv(Sw)*(u1-u2)';
% b) 画出投影前的两类样本点和投影线.(即特征向量：各类方差)
x1=m1(1:7,1);y1=m1(1:7,2);z1=m1(1:7,3);
x2=m2(1:7,1);y2=m2(1:7,2);z2=m2(1:7,3);
subplot(1,2,1);plot3(x1,y1,z1,'r*',x2,y2,z2,'b*'); %plot3为三维线图
title('原样本分布图');%画出原两类样本点
hold on
x=-2:0.1:2;
y=w(2,1)/w(1,1)*x;
z=w(3,1)/w(1,1)*x;
plot3(x,y,z);%画出投影线
legend('第一类样本','第二类样本','投影线',2);
%rotate3D on       %鼠标移动可以旋转三维坐标

% c) 画出投影后的两类样本点和投影线
subplot(1,2,2);
x=-2:0.1:2;
y=w(2,1)/w(1,1)*x;
z=w(3,1)/w(1,1)*x;
plot3(x,y,z);%画出投影线
hold on
% 对w1中的点投影
for i=1:7
    x1=(m1(i,1)+w(2,1)/w(1,1)*m1(i,2)+w(3,1)/w(1,1)*m1(i,3))/(1+(w(2,1)/w(1,1))^2+(w(3,1)/w(1,1))^2);
    y1=w(2,1)/w(1,1)*x1;
    z1=w(3,1)/w(1,1)*x1;
    plot3(x1,y1,z1,'r*');  
    X1(i,1)=x1;
    X1(i,2)=y1;
    X1(i,3)=z1;
end

% 对w2中的点投影
for i=1:7
    x2=(m2(i,1)+w(2,1)/w(1,1)*m2(i,2)+w(3,1)/w(1,1)*m2(i,3))/(1+(w(2,1)/w(1,1))^2+(w(3,1)/w(1,1))^2);
    y2=w(2,1)/w(1,1)*x2;
    z2=w(3,1)/w(1,1)*x2;
    plot3(x2,y2,z2,'b*');  title('投影后的样本分布图');
    X2(i,1)=x2;
    X2(i,2)=y2;
    X2(i,3)=z2;
end

for i=1:7
    X11(i)=w'*X1(i,:)';
    X21(i)=w'*X2(i,:)';
end
M1=mean(X11);
M2=mean(X21);
F=(M1+M2)/2;%分界阈值点
a=input('请输入测试样本点坐标  x=');
b=input('请输入测试样本点坐标  y=');
c=input('请输入测试样本点坐标  z=');
plot3(a,b,c,'go');
a1=(a+w(2,1)/w(1,1)*b+w(3,1)/w(1,1)*c)/(1+(w(2,1)/w(1,1))^2+(w(3,1)/w(1,1))^2);
b1=w(2,1)/w(1,1)*a1;
c1=w(3,1)/w(1,1)*a1;
L=[a1 b1 c1];
exam=w'*L';
if exam>F
   plot3(a1,b1,c1,'go'); 
   fprintf('该点属于第一类\n');
else
   plot3(a1,b1,c1,'go');
   fprintf('该点属于第二类\n');
end
rotate3D on    %鼠标移动可以旋转三维坐标