判别分析的MATLAB实现案例

%--------------------------------------------------------------------------%读取examp10_01.xls中数据，进行距离鉴别%--------------------------------------------------------------------------%********************************读取数据***********************************%读取文献examp10_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);%运用所创立的朴素贝叶斯分类器对象对训练样本进行鉴别，返回鉴别成果pre0，pre0也是字符串元胞向量pre0=ObjBayes.predict(meas);%运用confusionmat函数，并根据species和pre0创立混淆矩阵（包括总的分类信息的矩阵）[CLMat,order]=confusionmat(species,pre0);%以元胞数组形式查看混淆矩阵[[{'From/To'},order'];order,num2cell(CLMat)]%查看误判样品编号gindex1=grp2idx(pre0);%根据分组变量pre0生成一种索引向量gindex1gindex2=grp2idx(species);%根据分组变量species生成一种索引向量gindex2errid=find(gindex1~=gindex2)%通过对比两个索引向量，返回误判样品的观测序号向量%查看误判样品的误判状况head1={'Obj','From','To'};%设置表头%用num2cell函数将误判样品的观测序号向量errid转为元胞向量，然后以元胞数组形式查看误判成果[head1;num2cell(errid),species(errid),pre0(errid)]%对未知类别样品进行鉴别%定义未判样品观测值矩阵xx=[5.8 2.7 1.8 0.735.6 3.1 3.8 1.86.1 2.5 4.7 1.16.1 2.6 5.7 1.95.1 3.1 6.5 0.625.8 3.7 3.9 0.135.7 2.7 1.1 0.126.4 3.2 2.4 1.66.7 3 1.9 1.16.8 3.5 7.9 1];%运用所创立的朴素贝叶斯分类器对象对未判样品进行鉴别，返回鉴别成果pre1，pre1也是字符串元胞向量pre1=ObjBayes.predict(x)%--------------------------------------------------------------------------%加载fisheriris.mat中数据，进行Fisher鉴别%--------------------------------------------------------------------------%********************************加载数据***********************************loadfisheriris%把文献fisheriris.mat中数据导入MATLAB工作空间%**********************************待判样品*********************************%定义待判样品观测值矩阵xx=[5.8 2.7 1.8 0.735.6 3.1 3.8 1.86.1 2.5 4.7 1.16.1 2.6 5.7 1.95.1 3.1 6.5 0.625.8 3.7 3.9 0.135.7 2.7 1.1 0.126.4 3.2 2.4 1.66.7 3 1.9 1.16.8 3.5 7.9 1];%*********************************Fisher鉴别********************************%运用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);%提取各类的鉴别式得分ts1=ts(ts(:,1)==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类');%加标注框xlabel('第一鉴别式得分');%给X轴加标签ylabel('第二鉴别式得分');%给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构造鉴别式，%运用所有鉴别式看待判样品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(...)将所用鉴别式作用%在各组的组均值上，得到组均值投影矩阵，以表格形式返回这个矩阵.TabG是一种元胞%数组，形如%'Group''can1''can2'%'setosa'[-1.3849][1.8636]%'versicolor'[0.9892][1.6081]%'virginica'[1.9852][1.9443]%[class,TabCan,TabL,TabCon,TabM,TabG,trainscore]=fisher(...)返回%训练样品所对应的鉴别式得分trainscore.trainscore的第一列为各训练样品原本所%属类的类序号，第i+1列为第i个鉴别式得分.%%Copyrightxiezhh.%$Revision:1.0.0.0$$Date:/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.