%*************Fisher判別式**************
%*************FDA主元分析***************

%給定樣本集合
% X0=[5.1 3.5 0.4 0.2;4.9 3.0 1.4 0.2;4.7 3.2 1.3 0.2;4.6 3.1 1.5 0.2;5.0 3.6 1.4 0.2;5.4 3.9 1.7 0.4;4.6 3.4 1.4 0.3;5.0 3.4 1.5 0.2;4.4 2.9 1.4 0.2;4.9 3.1 1.5 0.1;5.4 3.7 1.5 0.2;4.8 3.4 1.6 0.2;4.8 3.0 1.1 0.1;4.3 3.0 1.1 0.1;5.8 4.0 1.2 0.2;5.7 4.4 1.5 0.4;5.4 3.9 1.3 0.4;5.1 3.5 1.4 0.3;5.7 3.8 1.8 0.3;5.1 3.8 1.5 5.3;5.4 3.4 1.7 0.2;5.1 3.7 1.5 0.4;4.6 3.6 1.0 0.2;5.1 3.3 1.7 0.5;4.8 3.4 1.9 0.2;5.0 3.0 1.6 0.2;5.0 3.4 1.4 0.2;5.2 3.5 1.5 0.2;5.2 3.4 1.4 0.2;4.7 3.2 1.6 0.2;4.8 3.1 1.6 0.2;5.4 3.4 1.5 0.4;5.2 4.1 1.5 0.1;5.5 4.2 1.4 0.2;4.9 3.1 1.5 0.2;5.0 3.2 1.2 0.2;5.5 3.5 1.3 0.2;4.9 3.6 1.4 0.1;4.4 3.0 1.3 0.2;5.1 3.4 1.5 0.2;5.0 3.5 1.3 0.3;4.5 2.3 1.3 0.3;4.4 3.2 1.3 0.2;5.0  3.5 1.6 0.6;5.1 3.8 1.9 0.3;4.8 3.0 1.4 0.3;5.1 3.8 1.6 0.2;4.6 3.2 1.4 0.2;5.3 3.7 1.5 0.2;5.0 3.3 1.4 0.3;7.0 3.2 4.7 1.4;6.4 3.2 4.5 1.5;6.9 3.1 4.9 1.5;5.5 2.3 4.0 1.3;6.5 2.8 4.6 1.5;5.7 2.8 4.5 1.3;6.3 3.3 4.7 1.6;4.9 2.4 3.3 1.0;5.6 2.9 4.6 1.3;5.2 2.7 3.9 1.4;5.0 2.0 3.5 1.0;5.9 3.0 4.2 1.5;6.0 2.2 4.0 1.0;6.1 2.9 4.7 1.4;5.6 2.9 3.9 1.3;6.7 3.1 4.4 1.4;5.6 3.0 4.5 1.5;5.8 2.7 4.1 1.0;6.2 2.2 4.5 1.5;5.6 2.5 3.9 1.1;5.9 3.2 4.8 1.8;6.1 2.8 4.0 1.3;6.3 2.5 4.9 1.5;6.1 2.8 4.7 1.2;6.4 2.9 4.3 1.3;6.6 3.0 4.4 1.4;6.8 2.8 4.8 1.4;6.7 3.0 5.0 1.7;6.0 2.9 4.5 1.5;5.7 2.6 3.5 1.0;5.5 2.4 3.8 1.1;5.5 2.4 3.7 1.0;5.8 2.7 3.9 1.2;6.0 2.7 5.1 1.6;5.4 3.0 4.5 1.5;6.0 3.4 4.5 1.6;6.7 3.1 4.7 1.5;6.3 2.3 4.4 1.3;5.6 3.0 4.1 1.3;5.5 2.5 5.0 1.3;5.5 2.6 4.4 1.2;5.7 3.0 4.6 1.4;5.8 2.6 4.0 1.2;5.0 2.3 3.3 1.0;5.6 2.7 4.2 1.3;5.7 3.0 4.2 1.2;5.7 2.9 4.2 1.3;6.2 2.9 4.3 1.3;5.1 2.5 3.0 1.1;5.7 2.8 4.1 1.3;6.3 3.3 6.0 2.5;5.8 2.7 5.1 1.9;7.1 3.0 5.9 2.1;5.3 2.9 5.6 1.8;6.5 3.0 5.8 2.2;7.6 3.0 5.6 2.1;4.9 2.5 4.5 1.7;7.3 2.9 6.3 1.8;6.7 2.5 5.8 1.8;7.2 3.6 6.1 2.5;6.5 3.2 5.1 2.0;6.4 2.7 5.3 1.9;6.8 3.0 5.5 2.1;5.7 2.5 5.0 2.0;5.8 2.8 5.1 2.4;5.4 3.2 5.3 2.3;6.5 3.0 5.5 1.8;7.7 3.8 6.7 2.2;7.7 2.6 5.9 2.3;5.0 2.2 5.0 1.5;5.9 3.2 5.7 2.3;5.6 2.8 4.9 2.0;7.7 2.8 6.7 2.0;6.3 2.7 4.9 1.8;6.7 3.3 5.7 2.1;7.2 3.2 6.0 1.8;6.2 2.8 4.8 1.8;6.1 3.0 4.9 1.8;6.4 2.8 5.6 2.1;7.2 3.0 5.8 1.6;7.4 2.8 6.1 1.9;7.9 3.8 6.4 2.0;6.4 2.8 5.6 2.2;6.3 2.8 5.1 1.5;6.1 2.6 5.6 1.4;7.7 3.0 6.1 2.3;6.3 3.4 5.6 2.4;6.4 3.1 5.5 1.8;6.0 3.0 4.8 1.8;6.9 3.1 5.4 2.1;6.7 3.1 5.6 2.4;6.9 3.1 5.1 2.3;5.8 2.7 5.1 1.9;6.8 3.2 5.9 2.3;6.7 3.3 5.7 2.5;6.7 3.0 5.2 2.3;6.3 2.5 5.0 1.9;6.5 3.0 5.2 2.0;6.2 3.4 5.4 2.3;5.9 3.0 5.1 1.8;]';
% X1=[1 1 1 2]';
% 
% X=[X0 X1 ];
% 
% %樣本元胞數(shù)組，元胞下標(biāo)與故障序號對應(yīng)
% XC{1,1}=X0;XC{1,2}=X1;   

run Iria;
XAve=(mean(X'))';

%初始化
Sb=zeros(m,m);
Sw=zeros(m,m);
Swk=zeros(m,m,c);

for j=1:c;
    Temp1=XC{1,j};
    k=length(Temp1(1,:));
    %樣本標(biāo)準(zhǔn)化XC
    for i=1:k;
        XC{1,j}(:,i)=(Temp1(:,i)-XAve)/std(Temp1(:,i));
    end
     P(j)=length(Temp1(1,:))/length(X(1,:));

    %類間離散度矩陣Sb
    XMean(:,j)=(mean(XC{1,j}'))';
    Sb=Sb+P(j)*(XMean(:,j)-XAve)*(XMean(:,j)-XAve)';
    %類內(nèi)離散度矩陣Sw
    for i=1:k;
        Temp2=XC{1,j}(:,i);
        Swk(:,:,j)=Swk(:,:,j)+(Temp2-XMean(:,j))*(Temp2-XMean(:,j))';
    end
    Sw=Sw+P(j)*Swk(:,:,j);
end
XCC=XC;

disp('類內(nèi)離散度矩陣')
Sb
disp('類內(nèi)離散度矩陣')
Sw

%總體散度矩陣St,及其特征值與特征向量
disp('總體散度矩陣')
% St=inv(Sw)*Sb
St=Sw+Sb;
[V,D]=eig(St);
disp('特征值和特征向量')
StD=diag(D)'%抽取相似對角陣D的對角向量
StV=V
disp('離散度準(zhǔn)則函數(shù)')
J=trace(St)

%特征值排序，提取主元
Stdmax=0;
for d=1:m-1;
    if StD(d+1)>StD(d)
        midD=StD(d+1);
        StD(d+1)=StD(d);
        StD(d)=midD;
        midV=StV(d+1);
        StV(d+1)=StV(d);
        StV(d)=midV;
        %樣本XC按主元排序
        for j=1:c;    
            Temp3=XC{1,j}(d+1,:);
            XC{1,j}(d+1,:)=XC{1,j}(d,:);
            XC{1,j}(d,:)=Temp3;
        end;
    end;
end;

disp('按序排列的特征值和特征向量')
StD
StV
disp('提取主元后的離散度準(zhǔn)則函數(shù)')
J=StD(1)+StD(2)

% for j=1:c;
%     disp(['第', int2str(j) , '分類主元空間'])
%     XC{1,j}(1:2,:)
% end

%第一、第二主元模式分布圖
plot(XC{1,1}(1,:),XC{1,1}(2,:),'b*')
title('基于FDA的Fisher數(shù)據(jù)分類')
hold on
plot(XC{1,2}(1,:),XC{1,2}(2,:),'rp')
hold on
plot(XC{1,3}(1,:),XC{1,3}(2,:),'k.')
hold on