fisher3d.m

能够实现Fisher线性判别

function y = fisher3d(w1,w2)

close all;clear;clc;
ts=cputime;
A = importdata('D:\MATLAB\R2007b\work\Pattern Recognition\fisher_testdata_3d.txt');
% 将实验数据存储为txt文件格式，用importdata命令读入数据
i=2;
dataset1=A( (i-1)*10+1:i*10 ,  4:6 );    % 将待分类的两类数据存入判别函数的输入参数x1,x2中
dataset2=A( (i-1)*10+1:i*10 ,  7:9 );
dataset3=A( (i-1)*10+1:i*10 ,  1:3 );

% dataset1=randn(10,3)-2*[ones(1,10);ones(1,5),zeros(1,5);zeros(1,5),ones(1,5)]';  % 利用随机数据来检测程序的分类性能
% dataset2=1*[ones(1,10);ones(1,5),zeros(1,5);zeros(1,5),ones(1,5)]'+randn(10,3);

[n1,d1]=size(dataset1);   % 读取输入数据的个数n和维数d
[n2,d2]=size(dataset2);

% dataset1w=[dataset1,1*ones(1,n1)'];     % 对数据作出类别标记
% dataset2w=[dataset1,2*ones(1,n2)'];
% dataset3w=[dataset1,3*ones(1,n3)'];

if d1~=d2           % 如果维数不等，则提示错误
    error('The dimensions of w1 and w2 must be the same !');
    return;
end
if d1>3
    error('The dimensions of the data is too large to be shown !');
    return;
end

mw1=mean(dataset1);     % 求两类数据的类内均值
mw2=mean(dataset2);

s1=0;   s2=0;       % 求两类数据的类内散布矩阵si和总类内散布矩阵sw
for i=1:n1
    s1=s1+(dataset1(i,:)'-mw1')*(dataset1(i,:)'-mw1')';
end
for i=1:n2
    s2=s2+(dataset2(i,:)'-mw2')*(dataset2(i,:)'-mw2')';
end
sw=s1+s2;

w=sw^(-1)*(mw1'-mw2')   % 求出最优方向w

yw1=w'*dataset1';          % y1是第一类数据在最优方向上的投影点的值
yw2=w'*dataset2';          % y2是第二类数据在最优方向上的投影点的值
yw=[yw1,yw2];         % y0是两类数据在最优方向上的投影点的值
% yalpha=mean(y0);    % 判别阈值设为点集{y0}的平均值
% y=y0-yalpha;        % y是经过判别处理后数据的投影值，y>0则判别为第一类，反之判为第二类

x1=dataset1(:,1);y1=dataset1(:,2);z1=dataset1(:,3);     % 将各点数据分为x,y,z三维数据，方便作图
x2=dataset2(:,1);y2=dataset2(:,2);z2=dataset2(:,3);
x=[x1;x2];y=[y1;y2];z=[z1;z2];
% xmax=max(abs(x));ymax=max(abs(y));zmax=max(abs(z));
% axismax=max([xmax;ymax;zmax]);
axismax=max(abs(x));            % 求出数据点在x轴上的最大值
Aw=sqrt(w(1)^2+w(2)^2+w(3)^2);  % 求出最优方向W的幅值
ye=1.5*axismax/Aw;      % ye,ys分别是x轴的刻度范围
ys=-ye;

figure;
subplot(121);
    stem(yw);Grid;  % 画出分类结果

% figure;
subplot(122);             % 在二维平面上画出分类结果
    h1=stem3(x1,y1,z1,'bo');            % 画出第一类数据点
        set(h1,'LineStyle','none');     
        set(h1,'MarkerFace','red');
        hold on;
    h2=stem3(x2,y2,z2,'bs');            % 画出第二类数据点
        set(h2,'LineStyle','none');     
        set(h2,'MarkerFace','green');
        hold on;
    h3=stem3(ye*w(1),ye*w(2),ye*w(3),'md');     % 画出最优方向W直线的正方向箭头
        set(h3,'LineStyle','none');  
        set(h3,'MarkerFace','blue');
        hold on;
    h4=line([0 ye*w(1)],[0 ye*w(2)],[0 ye*w(3)],'LineWidth',2,'color',[1 0 1]);     % 画出最优方向W直线，正半轴用粗线表示
    h5=line([ys*w(1) 0],[ys*w(2) 0],[ys*w(3) 0],'LineWidth',1,'color',[0 0 1]);     % 负半轴用细线表示
    h6=line([0 0],[0 0],[0 -1],'LineWidth',2,'color',[0 1 0]);                      % 画出通过原点与Z轴平行的参考线

pjpoint1=[];        % 画出各点到最优方向直线的投影线，计算出投影点
for i=1:n1
    p1=projectline3d([x1(i) y1(i) z1(i)],w,'r');
    pjpoint1=[pjpoint1,p1];
end
pjpoint2=[];
for j=1:n2
    p2=projectline3d([x2(j) y2(j) z2(j)],w,'g');
    pjpoint2=[pjpoint2,p2];
end
pjpoint=[pjpoint1,pjpoint2];  

xlabel('x(1)--axis x');ylabel('x(2)--axis y');zlabel('x(3)--axis z');   % 标出坐标轴
axis equal;
te=cputime-ts
cost=whos;                 %将内存信息读入cost 
sumcost=0; 
for i=1:length(cost)       %一下子没找到不用循环的方法，先凑合着用 
   sumcost=sumcost+cost(i).bytes; 
end 
clear cost;                %cost变量本身占用一定空间，此时清除它 
sumcost