MATLAB数据分析方法-(3)ppt课件

普通高等院校计算机课程规划教材,MATLAB数据分析方法,李柏年吴礼斌主编张孔生丁华参编,第4章判别分析,判别分析的基本思想是根据已知类别的样本所提供的信息，总结出分类的规律性，建立判别公式和判别准则，判别新的样本点所属类型。本章介绍距离判别分析、Bayes判别分析极其MATLAB软件的实现。,4.1距离判别分析,4.1.1判别分析的概念,在一些自然科学和社会科学的研究中，研究对象用某种方法已划分为若干类型，当得到的一个新样品数据（通常是多元的），要确定该样品属于已知类型中的哪一类，这样的问题属于判别分析.,从统计数据分析的角度，可概括为如下模型：,设有k个总体,它们都是p元总体,其数量指标是1)若总体的分布函数是已知，对于任一新样品数据,判断它来自哪一个总体。,2)通常各个总体的分布是未知的，由从各个总体取得的样本（训练样本）来估计。一般，先估计各个总体的均值向量与协方差矩阵。,称称为n维向量x，y之间的闵可夫斯基距离，其中为常数。,欧氏距离,显然，当r=2和1时闵可夫斯基距离分别为欧氏距离和绝对距离.,(3)两个总体之间的马氏距离,设有两个总体G1,G2，两个总体的均值向量分别为，协方差矩阵相等，皆为,则两个总体之间的马氏距离为,通常，在判别分析时不采用欧氏距离的原因在于，该距离与量纲有关.例如平面上有A,B,C,D四个点，横坐标为代表重量（单位：kg），纵坐标代表长度（单位：cm），如下页图。,(4.1.3),这时,显然ABCD,如果现在长度用mm为单位，重量的单位保持不变，于是A点的坐标为(0,50)，B点的坐标为(0,100)，此时计算线段的长度为,此时，AB0，x属于G1;若W(x)0);p0=(n11+n22)/(n1+n2)%计算回代误判率,fori=1:n1A=a(1:i-1,i+1:n1,:);n1=length(A(:,1);n2=length(b(:,1);s1=cov(A);s2=cov(b);p=4;s=(n1-1)*s1+(n2-1)*s2)/(n1+n2-2);D11(i)=(a(i,:)-mean(A)*inv(s)*(a(i,:)-mean(A)-(a(i,:)-mean(b)*inv(s)*(a(i,:)-mean(b);endfori=1:n2B=b(1:i-1,i+1:n2,:);n1=length(a(:,1);n2=length(B(:,1);s1=cov(A);s2=cov(B);p=4;s=(n1-1)*s1+(n2-1)*s2)/(n1+n2-2);D22(i)=(b(i,:)-mean(B)*inv(s)*(b(i,:)-mean(B)-(b(i,:)-mean(a)*inv(s)*(b(i,:)-mean(a);endN11=length(find(D110);N22=length(find(D220);p1=(N11+N22)/(n1+n2)%计算交叉误判率,输出结果：p0=0.1923p1=0.2400,4.3Bayes判别分析,贝叶斯公式是一个我们熟知的公式,距离判别只要求知道总体的数字特征，不涉及总体的分布函数，当参数和协方差未知时，就用样本的均值和协方差矩阵来估计。距离判别方法简单实用，但没有考虑到每个总体出现的机会大小，即先验概率，没有考虑到错判的损失。贝叶斯判别法正是为了解决这两个问题提出的判别分析方法。,因此，两个总体的Bayes判别准则为,2.两个正态总体的Bayes判别,（1）两个总体协方差矩阵相等的情形,设总体G1,G2的协方差矩阵相等且为，概率密度函数为：,(4.3.2),损失相等的Bayes判别准则为,其中,基于两正态总体后验概率的Bayes判别准则为,其中,在实际问题中，关于先验概率，通常用下列两种方式选取:1）采用等概率选取，即2）按训练样本的容量的比例选取，即,例4.3.1对例4.1.1的数据，重新对上述三个蠓虫的类别进行Bayes判别.（假设误判损失相等）解：第1步：可以验证两个总体服从二元正态分布；（第二章的正态性检验，读者自证）第2步：检验两个总体的协方差矩阵相等；第3步：估计两个总体的先验概率,这里按样本容量的比例选取.由于Apf与Af分别为6个与9个，故估计Apf类蠓虫的先验概率，Af类蠓虫的先验概率；,第4步：利用MATLAB软件计算：,apf=1.14,1.78;1.18,1.96;1.20,1.86;1.26,2.;1.28,2;1.30,1.96;af=1.24,1.72;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90;1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08;x=1.24,1.8;1.28,1.84;1.4,2.04;m1=mean(apf);m2=mean(af);s1=cov(apf);s2=cov(af);s=(5*s1+8*s2)/13;fori=1:3w1(i)=m1*inv(s)*x(i,:)-1/2*m1*inv(s)*m1+log(0.4);w2(i)=m2*inv(s)*x(i,:)-1/2*m2*inv(s)*m2+log(0.6);ifw1(i)=w2(i)disp(第,num2str(i),个蠓虫属于Apf类);elsedisp(第,num2str(i),个蠓虫属于Af类);end;end;输出结果：第1个蠓虫属于Apf类第2个蠓虫属于Apf类第3个蠓虫属于Apf类,（2）两个总体协方差矩阵不相等的情形,设总体的协方差矩阵不相等分别为1,2概率密度函数为：,则基于两正态总体误判损失相等的Bayes判别准则,其中,例4.3.2对破产的企业收集它们在破产前两年的年度财务数据，对财务良好的企业也收集同一时间的数据.数据涉及四个变量:现金流量/总债务，净收益/总资产，流动资产/流动债务，以及流动资产/净销售额，数据如表4.2所示.假定两总体G1,G2均服从四元正态分布，在误判损失相等且先验概率按比例分配的条件下，对待判样本进行bayes判别.表4.2,解：第1步：检验两个总体的协方差矩阵相等；源程序如下：,A=-0.45-0.411.090.450.510.102.490.54-0.13-0.141.420.440.170.071.800.52x=-0.23-0.300.330.18;0.150.052.170.55-0.28-0.231.190.66;0.480.091.240.18;G1=A(:,1:4);G2=A(:,5:8);%二类总体数据m1=mean(G1);m2=mean(G2);s1=cov(G1);s2=cov(G2);n=18;n1=9;n2=9;p=2;s=(n1-1)*s1+(n2-1)*s2)/(n1+n2-2);Q1=(n1-1)*(log(det(s)-log(det(s1)-p+trace(inv(s)*s1);Q2=(n2-1)*(log(det(s)-log(det(s2)-p+trace(inv(s)*s2);ifQ1chi2inv(0.95,p*(p+1)/2)输出结果：两组数据协方差不全相等,第2步：根据第1步结论，构造判别函数，得出判结果.p1=n1/n;p2=n2/n;%计算先验概率fori=1:4d1(i)=mahal(x(i,:),G1)-log(det(s1)-2*log(p1);d2(i)=mahal(x(i,:),G2)-log(det(s2)-2*log(p2);ifd1(i)=d2(i)disp(第,num2str(i),个属于破产企业);elsedisp(第,num2str(i),个属于非破产企业);end;end;输出结果：第1个属于破产企业第2个属于非破产企业第3个属于破产企业第4个属于非破产企业,4.3.2多个总体的Bayes判别,设有k个总体G1,G2,Gk的概率密度为fj(x)各总体出现的先验概率为,1.一般讨论,当出现样品,时,总体,的后验概率,Bayes判别准则为:若,则判样本,注：当达到最大后验概率的,不止一个时，可判,为达到最大后验概率的总体的任何一个.,2.多个正态总体的Bayes判别,（1）当,时，设,线性判别函数为,其中,基于误判损失相等的Bayes判别准则为,基于后验概率的Bayes判别准则为,其中,在实际问题中，由于未知，各总体的训练样本均值（2）当,不全相等时,设,则基于后验概率的Bayes判别准则为,其中,未知，,估计.,例4.3.3.某医院利用心电图检测来对人群进行划分，数据见表.“g=1”表示健康人，“g=2”表示主动脉硬化患者，“g=3”表示冠心病患者，X1,X2表示测得的心电图中表明心脏功能的两项不相关的指标.某受试者心电图该两项指标的数据分别为380.20，9.08.设先验概率按比例分配，进行bayes判别，判定其归属.,表4.324人心电图数据,解：A=261.017.36185.395.99189.595.46x=380.209.08;G1=A(1:11,:);G2=A(12:18,:);G3=A(19:23,:);%三类总体数据n=23;k=3;p=2;n1=11;n2=7;n3=5;f=p*(p+1)*(k-1)/2;d=(2*p2+3*p-1)*(1/(n1-1)+1/(n2-1)+1/(n3-1)-1/(n-k)/(6*(p+1)*(k-1);p1=n1/n;p2=n2/n;p3=n3/n;m1=mean(G1);m2=mean(G2);m3=mean(G3);s1=cov(G1);s2=cov(G2);s3=cov(G3);%计算协方差阵s=(n1-1)*s1+(n2-1)*s2+(n3-1)*s3)/(n-k);M=(n-k)*log(det(s)-(n1-1)*log(det(s1)+(n2-1)*log(det(s2)+(n3-1)*log(det(s3);T=(1-d)*M%计算统计量观测值,C=chi2inv(0.95,f)ifTchi2inv(0.95,f)disp(三组数据协方差相等);elsedisp(三组数据协方差不全相等);end;w(1)=m1*inv(s)*x-1/2*m1*inv(s)*m1+log(p1);w(2)=m2*inv(s)*x-1/2*m2*inv(s)*m2+log(p2);w(3)=m3*inv(s)*x-1/2*m3*inv(s)*m3+log(p3);fori=1:3ifw(i)=max(w)disp(属于第,num2str(i),组);end;end;输出结果:三组数据协方差相等属于第2组,4.3.3平均误判率,Byaes判别的有效性可以通过平均误判率来确定。这里仅对两个正态总体,且协方差矩阵相等,的情况下研究平均误判率的计算.,设总体,，其先验概率,，两个总体,的,马氏平方距离记为,则基于误判损失相等时的平均误判率为,其中,为标准正态分布函数.,从（4.3.11）式知，当总体,的马氏平方距离,越大，即两总体的分离程度越大时，平均误判概率最小.推广到一般情况也成立.,（4.3.11）,例4.3.42008年全国部分地区城镇居民人均年家收入情况见表.按四种指标分为二类，用bayes判别判定青海、广东两省区属于哪一类，并用回代法和交叉法对误判率进行估计(假定误判损失相等).表4.4,解：第1步，检验三个总体的协方差矩阵相等；A=18738.96778.36452.757707.879422.22938.15141.751976.49;x=8595.48763.0750.173458.6315188.392405.92701.253382.95;%待判样品G1=A(1:2,:);G2=A(3:8,:);G3=A(9:27,:);%输入三类总体数据n1=size(G1,1);%总体G1的样本数n2=size(G2,1);%总体G2的样本数n3=size(G3,1);%总体G3的样本数,n=n1+n2+n3;%三个总体合并的样本数k=3;p=4;f=p*(p+1)*(k-1)/2;d=(2*p2+3*p-1)*(1/(n1-1)+1/(n2-1)+1/(n3-1)-1/(n-k)/(6*(p+1)*(k-1);p1=n1/n;p2=n2/n;p3=n3/n;m1=mean(G1);m2=mean(G2);m3=mean(G3);s1=cov(G1);s2=cov(G2);s3=cov(G3);%计算协方差阵s=(n1-1)*s1+(n2-1)*s2+(n3-1)*s3)/(n-k);M=(n-k)*log(det(s)-(n1-1)*log(det(s1)+(n2-1)*log(det(s2)+(n3-1)*log(det(s3);T=(1-d)*M%计算统计量观测值C=chi2inv(0.95,f)ifTchi2inv(0.95,f)disp(三组数据协方差相等);elsedisp(三组数据协方差不全相等);end,输出结果：三组数据协方差相等,第2步，根据第1步结论，构造判别函数，得出判别结果.fori=1:2w(1)=m1*inv(s)*x(i,:)-1/2*m1*inv(s)*m1+log(p1);w(2)=m2*inv(s)*x(i,:)-1/2*m2*inv(s)*m2+log(p2);w(3)=m3*inv(s)*x(i,:)-1/2*m3*inv(s)*m3+log(p3);%计算判别函数forj=1:3ifw(j)=max(w)disp(待判样品属于第,num2str(j),类城市);endendend输出结果：待判样品属于第3类城市待判样品属于第2类城市,第3步，计算回代误判率.n11=0;n22=0;n33=0;,fori=1:n1w1(i,1)=m1*inv(s)*G1(i,:)-1/2*m1*inv(s)*m1+log(p1);w1(i,2)=m2*inv(s)*G1(i,:)-1/2*m2*inv(s)*m2+log(p2);w1(i,3)=m3*inv(s)*G1(i,:)-1/2*m3*inv(s)*m3+log(p3);%计算判别函数forj=1:3ifw1(i,j)=max(w1(i,:)endendend,fori=1:n3w3(i,1)=m1*inv(s)*G3(i,:)-1/2*m1*inv(s)*m1+log(p1);w3(i,2)=m2*inv(s)*G3(i,:)-1/2*m2*inv(s)*m2+log(p2);w3(i,3)=m3*inv(s)*G3(i,:)-1/2*m3*inv(s)*m3+log(p3);%计算判别函数forj=1:3ifw3(i,j)=max(w3(i,:)endendendp00=(n11+n22+n33)/(n1+n2+n3),输出结果：p00=0,第4步，计算交叉误判率.,N11=0;N22=0;N33=0;,fork=1:n1A=G1(1:k-1,k+1:n1,:);N1=length(A(:,1);M1=mean(A,1);s11=cov(A);S1=(N1-1)*s11+(n2-1)*s2+(n3-1)*s3)/(N1+n2+n3-k);P01=N1/(n-1);P02=n2/(n-1);P03=n3/(n-1);%计算先验概率fori=1:n1W1(i,1)=M1*inv(S1)*G1(i,:)-1/2*M1*inv(S1)*M1+log(P01);W1(i,2)=m2*inv(S1)*G1(i,:)-1/2*m2*inv(S1)*m2+log(P02);W1(i,3)=m3*inv(S1)*G1(i,:)-1/2*m3*inv(S1)*m3+log(P03);%计算判别函数forj=1:3ifW1(i,j)=max(W1(i,:)endendendend,fork=1:n2B=G2(1:k-1,k+1:n2,:);N2=length(B(:,1);M2=mean(B,1);s22=cov(B);S2=(n1-1)*s1+(N2-1)*s22+(n3-1)*s3)/(n1+N2+n3-k);%计算混合