matlab之判别分析

      想象一下，你是一个土豪，有很多企业来找你融资，你不想投资给将来会破产的企业（那多亏啊！还不如给女票呢），可是，你又不知道哪个企业会破产，于是，你找到了gbl来帮你判断，gbl听说你的困难之后，即兴讲解了一番，解决了你的难题！

     这类问题属于判别分析，首先，找到一些运行良好的企业，搜集他的指标，取为类G1，然后找到破产的企业，搜集他的指标，取为类G2，最后，就是来判断找你融资的企业x属于G1还是G2？属于G1就赶紧掏钱啊！

     考虑一种最简单的情况，G1和G2的协方差阵相等。分别用距离判别，fisher判别，bayes判别来建立线性判别函数进行判别。

     在判别之前，先看一下判别分析有木有意义！！！

%两总体G1,G2,x1是G1的样本,x是合并的样本(判别效果分析)
x1=x(1:21,:);
x2=x(22:46,:);
m=4;%四项运行指标
n1=21;
n2=25;
bar1=mean(x1)';bar2=mean(x2)';%列向量
A1=cov(x1)*(n1-1);
A2=cov(x2)*(n2-1);
S=(A1+A2)/(n1+n2-2);%合并的样本协方差阵，并不等于cov(x)
d2=(bar1-bar2)'*inv(S)*(bar1-bar2);
f=(n1+n2-m-1)*n1*n2*d2/(m*(n1+n2)*(n1+n2-2));
p=1-fcdf(f,m,n1+n2-m-1);
if p<0.05
   disp('拒绝原假设,有显著差异，判别分析有意义')
else
    disp('接受原假设，无差异,判别分析没有意义')
end 

有意义的话，就可以进行判别分析咯！没有意义就没必要判别啦！

%距离判别线性判别函数,x中的47到54行是待判样品，1-21是第一类训练样品，22-46是第二类训练样品
clear;clc;
load('x.mat')
load('num.mat')%1-46号样品的类别，破产或运行良好
y1_x=(inv(S)*bar1)'%求系数
y1_c=bar1'*inv(S)*bar1
y2_x=(inv(S)*bar2)'
y2_c=bar2'*inv(S)*bar2
ans=[];
for i=1:54
   x0=x(i,:)';
   Y1=(inv(S)*bar1)'*x0-(1/2)*bar1'*inv(S)*bar1;
   Y2=(inv(S)*bar2)'*x0-(1/2)*bar2'*inv(S)*bar2;
   if Y1>Y2
       ans(i)=1;
   else
       ans(i)=2;
   end
end
huipan=ans(1:46)';%回判
sure=sum(ans(1:46)==num(1:46)')/46 %回判正确的比率
daipan=ans(47:54)’%待判样品类别

%fisher判别 
clear;clc;
load('x.mat')
load('num.mat')
ans=[];
A=A1+A2;
d2=(bar1-bar2)'*inv(A)*(bar1-bar2);
d=sqrt(d2);
l=inv(A)*(bar1-bar2)/d;
ubar=l'*(bar1+bar2)/2 %判别 
ux=inv(A)*(bar1-bar2)/d %线性判别函数
for i=1:54
   x0=x(i,:)';
   u=x0'*inv(A)*(bar1-bar2)/d;
   if u>ubar
       ans(i)=1;
   else
       ans(i)=2;
   end
end
huipan=ans(1:46)';%回判
sure=sum(ans(1:46)==num(1:46)')/46
daipan=ans(47:54)'%待判

%Bayes线性判别函数，假设q1=q2=q

clear;clc;
load('x.mat')
load('num.mat')
ans=[];
q=0.5
y1_x=(inv(S)*bar1)'%求系数
y1_c=log(q)-(1/2)*bar1'*inv(S)*bar1
y2_x=(inv(S)*bar2)'
y2_c=log(q)-(1/2)*bar2'*inv(S)*bar2
for i=1:54
   x0=x(i,:)';
   %取大的
   Y1=(inv(S)*bar1)'*x0+log(q1)-(1/2)*bar1'*inv(S)*bar1;
   Y2=(inv(S)*bar2)'*x0+log(q2)-(1/2)*bar2'*inv(S)*bar2;
   if Y1>Y2
       ans(i)=1;
   else
        ans(i)=2;
   end
end
huipan=ans(1:46)';
sure=sum(ans(1:46)==num(1:46)')/46%回判正确率
daipan=ans(47:54)'