判别分析的MATLAB实现案例

% % 读取exampl0_01.xls中数据，进行距离判别% 、‘、丁-%读取文件exampl0_01.xls的第1个工作表中C2:F51范围的数据，即全部样本数据，包括未判企业sample=xlsread('examp10_01.xls','','C2:F51');%读取文件examp10_01.xls的第1个工作表中C2:F47范围的数据，即已知组别的样本数据,training=xlsread('examp10_01.xls','','C2:F47');%读取文件examp10_01.xls的第1个工作表中B2:B47范围的数据，即样本的分组信息数据,group=xlsread('examp10_01.xls','','B2:B47');obs=[1:50]';%企业的编号屮、%距离判别，判别函数类型为mahalanobis，返回判别结果向量C和误判概率err[C,err]=classify(sample,training,group,'mahalanobis');[obs,C]%查看判别结果err%查看误判概率% % 加载fisheriris.mat中数据，进行贝叶斯判别% ”£八口loadfisheriris%把文件fisheriris.mat中数据导入MATLAB工作空间head0={'Obj','x1','x2','x3','x4','Class'};%设置表头[head0;num2cell([[1:150]',meas]),species]%以元胞数组形式查看数据/介*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^%用meas和species作为训练样本，创建一个朴素贝叶斯分类器对象ObjBayesObjBayes=NaiveBayes.fit(meas,species);%利用所创建的朴素贝叶斯分类器对象对训练样本进行判别，返回判别结果preO,pre0也是字符串元胞向量pre0=ObjBayes.predict(meas);%利用confusionmat函数，并根据species和pre0创建混淆矩阵(包含总的分类信息的矩阵)[CLMat,order]=confusionmat(species,pre0);%以元胞数组形式查看混淆矩阵[[{'From/To'},order'];order,num2cell(CLMat)]%查看误判样品编号gindexl=grp2idx(pre0);%根据分组变量preO生成一个索引向量gindexlgindex2=grp2idx(species);%根据分组变量species生成一个索引向量gindex2errid=find(gindex1~=gindex2)%通过对比两个索引向量，返回误判样品的观测序号向量%查看误判样品的误判情况head1={'Obj','From','To'};%设置表头%用num2cell函数将误判样品的观测序号向量errid转为元胞向量，然后以元胞数组形式查看误判结果[head1;num2cell(errid),species(errid),preO(errid)][5.82.71.8O.735.63.13.81.86.12.54.71.16.12.65.71.95.13.16.5O.625.83.73.9O.135.72.71.1O.126.43.22.41.66.731.91.16.83.57.91x对未知类别样品进行判别定义未判样品观测值矩阵x];利用所创建的朴素贝叶斯分类器对象对未判样品进行判别，返回判别结果prel,pre1也字符串元胞向量pre1=ObjBayes.predict(x)% % 加载fisheriris.mat中数据，进行Fisher判别% U/ 11 \匸floadfisheriris%把文件fisheriris.mat中数据导入MATLAB工作空间U/ _i-r*•~r*11A,—i-*>—i*4^*4^*4^*4^*4^*4^*4^*4^*4^*4^*4^*4^*4^*4^*4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>**4>*/AXj>XT»XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>XT>*TvZ|\| /*II/]—1—.t-1|—t%定义待判样品观测值矩阵xx=[5.82.71.80.735.63.13.81.86.12.54.71.16.12.65.71.95.13.16.50.625.83.73.90.135.72.71.10.126.43.22.41.66.731.91.16.83.57.91];以/n^f**^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^*^|iI iI{_ir—J—*iI—W|I%利用fisher函数进行判别，返回各种结果(见fisher函数的注释)[outclass,TabCan,TabL,TabCon,TabM,TabG]=fisher(x,meas,species)%************************绘制两个判别式得分的散点图**************************%利用fisher函数进行判别，返回各种结果，其中ts为判别式得分[outclass,TabCan,TabL,TabCon,TabM,TabG,ts]=fisher(x,meas,species);%提取各类的判别式得分tsl = ts(ts(:,l) == 1,:); % setosa类的判别式得分ts2 = ts(ts(:,1) == 2,:); % versicolor类的判别式得分ts3 = ts(ts(:,1) == 3,:); % virginica类的判别式得分plot(ts1(:,2),ts1(:,3),'ko') %setosa类的判别式得分的散点图holdonplot(ts2(:,2),ts2(:,3),'k*') %versicolor类的判别式得分的散点图plot(ts3(:,2),ts3(:,3),'kp') %virginica类的判别式得分的散点图legend('setosa类'，’versicolor类'，’virginica类')；％加标注框xlabelC第一判别式得分')；％给X轴加标签ylabelC第二判别式得分')；％给Y轴加标签%************************只用一个判别式进行Fisher判别************************%令fisher函数的第4个输入为0.5，就可以只用一个判别式进行判别[outclass,TabCan,TabL,TabCon,TabM,TabG]=fisher(x,meas,species,0.5)function [outclass,TabCan,TabL,TabCon,TabM,TabG,trainscore]fisher(sampledata,training,group,contri)%FISHER判别分析.%class=fisher(sampledata,training,group)根据训练样本training构造判另U式,% 利用所有判别式对待判样品sampledata进行判别.sampledata和training是具有相同% 列数的矩阵，它们的每一行对应一个观测，每一列对应一个变量group是training对% 应的分组变量，它的每一个元素定义了training中相应观测所属的类.group可以是% 个分类变量，数值向量，字符串数组或字符串元胞数组.training和group必须具有相%同的行数.fisher函数把group中的NaN或空字符串作为缺失数据，从而忽略training% 中相应的观测.class中的每个元素指定了sampledata中的相应观测所判归的类，它和% group具有相同的数据类型.%%class=fisher(sampledata,training,group,contri)根据累积贡献率不低于% contri,确定需要使用的判别式个数，默认情况下，使用所有判别式进行判别contri%是一个在(0,1]区间内取值的标量，用来指定累积贡献率的下限.%%[class,TabCan]=fisher(...)以表格形式返回所用判别式的系数向量，若contri%取值为1,则返回所有判别式的系数向量.TabCan是一个元胞数组，形如%'Variable''can1''can2'%'x1'[-0.2087][0.0065]%'x2'[-0.3862][0.5866]%'x3'[0.5540][-0.2526]%'x4'[0.7074][0.7695]%[class,TabCan,TabL]=fisher(...)以表格形式返回所有特征值，贡献率，累积%贡献率等.TabL是-个元胞数组，形如%'Eigenvalue''Difference''Proportion''Cumulative'%[32.1919][31.9065][0.9912][0.9912]%[0.2854][][0.0088][1]%%[class,TabCan,TabL,TabCon]=fisher(...)以表格形式返回混淆矩阵(包含总%的分类信息的矩阵).TabCon是个元胞数组，形如%'From/To''setosa''versicolor''virginica'%'setosa'[50][0][0]%'versicolor'[0][48][2]%'virginica'[0][1][49]%%[class,TabCan,TabL,TabCon,TabM]=fisher(...)以表格形式返回误判矩阵.%TabM是一个元胞数组，形如%'Obj''From''To'%[71]'versicolor''virginica'%[84]'versicolor''virginica'%[134]'virginica''versicolor'%

%[class,TabCan,TabL,TabCon,TabM,TabG]=fisher(...)将所用判别式作用%数组，形如%'Group''can1'%'setosa'[-1.3849]%'versicolor'[0.9892]%'virginica'[1.9852]%[class,TabCan,TabL,TabCon,%在各组的组均值上，得到组均值投影矩阵，以表格形式返回这个矩阵TabG%[class,TabCan,TabL,TabCon,'can2'[1.8636][1.6081][1.9443]TabM,TabG,trainscore]=fisher(...)返回% 训练样品所对应的判别式得分trainscore.trainscore的第一列为各训练样品原本所% 属类的类序号，第i+1列为第i个判别式得分.%Copyright2009xiezhh.%$Revision:1.0.0.0$$Date:2009/10/0310:40:34$ifnargin<3error（'错误：输入参数太少，至少需要3个输入•’）;end%根据分组变量生成索引向量gindex，组名元胞向量groups，组水平向量glevels[gindex,groups,glevels]=grp2idx(group);%忽略缺失数据nans=find(isnan(gindex));if~isempty(nans)training(nans,:)=[];gindex(nans)=[];endngroups=length(groups);gsize=hist(gindex,1:ngroups);nonemptygroups=find(gsize>0);nusedgroups=length(nonemptygroups);%判断是否有空的组ifngroups>nusedgroupswarning('警告：有空的组•')；end[n,d]=size(training);ifsize(gindex,1)~=nerror('错误：输入参数大小不匹配，GROUP与TRAINING必须具有相同的行数•’);elseifisempty(sampledata)sampledata=zeros(0,d,class(sampledata));elseifsize(sampledata,2)~=derror('错误：输入参数大小不匹配,SAMPLEDATA与TRAINING必须具有相同的列数end%设置contri的默认值为1,并限定contri在(0,1]内取值ifnargin<4||isempty(contri)contri=1;endif~isscalar(contri)||contri>1||contri<=0error('错误：contri必须是一个在(0,1]内取值的标量.’)；endifany(gsize==1)error('错误：TRAINING中的每个组至少应有两个观测.’)；end%计算各组的组均值gmeans=NaN(ngroups,d);fork=nonemptygroupsgmeans(k,：)=mean(training(gindex==k,：),1);end%计算总均值totalmean=mean(training,1);%计算组内离差平方和矩阵E和组间离差平方和矩阵BE=zeros(d);B=E;fork=nonemptygroups%分别估计各组的组内离差平方和矩阵.[Q,Rk]=qr(bsxfun(@minus,training(gindex==k,：),gmeans(k,：)),0);%各组的组内离差平方和矩阵：AkHat=Rk'*Rk%判断各组的组内离差平方和矩阵的正定性s=svd(Rk);ifany(s<=max(gsize(k),d)*eps(max(s)))error('错误：TRAINING中各组的组内离差平方和矩阵必须是正定矩阵');endE=E+Rk'*Rk;%计算总的组内离差平方和矩阵E%计算组间离差平方和矩阵BB=B+(gmeans(k,：)-totalmean)'*(gmeans(k,：)-totalmean)*gsize(k);end%求inv(E)*B的正特征值与相应的特征向量EB=E\B;[V,D]=eig(EB);D=diag(D);[D,idD]=sort(D,'descend');%将特征值按降序排列V=V(:,idD);NumPosi=min(ngroups-1,d);%确定正特征值个数D=D(1:NumPosi,:);CumCont=cumsum(D/sum(D));%计算累积贡献率%以表格形式返回所有特征值，贡献率，累积贡献率等.TabL是一个元胞数组head={'Eigenvalue','Difference','Proportion','Cumulative'};TabL=cell(NumPosi+1,4);TabL(1,:)=head;TabL(2:end,1)=num2cell(D);ifNumPosi==1TabL(2:end-1,2)={0};elseTabL(2:end-1,2)=num2cell(-diff(D));endTabL(2:end,3)=num2cell(D/sum(D));TabL(2:end,4)=num2cell(CumCont);%根据累积贡献率的下限contri确定需要使用的判别式个数CumContGeConCumContGeCon=find(CumCont>=contri);CumContGeCon=CumContGeCon(1);V=V(:,1:CumContGeCon);%需要使用的判别式系数矩阵%以表格形式返回所用判别式的系数向量，若contri取值为1,%则返回所有判别式的系数向量.TabCan是一个元胞数组TabCan=cell(d+1,CumContGeCon+1);TabCan(1,1)={'Variable'};TabCan(2:end,1)=strcat('x',cellstr(num2str((1:d)')));TabCan(1,2:end)=strcat('can',cellstr(num2str((1:CumContGeCon)')));TabCan(2:end,2:end)=num2cell(V);%将训练样品与待判样品放在一起进行判别m=size(sampledata,1);gv=gmeans*V;stv=[sampledata;training]*V;nstv=size(stv,1);message='';outclass=NaN(nstv,1);fori=1:nstvobji=bsxfun(@minus,stv(i,:),gv);obji=sum(obji."2,2);idclass=find(obji==min(obji));iflength(idclass)>1idclass=idclass(1);message='警告:出现