MATLAB数据分析方法

一般高等院校计算机课程规划教材MATLAB数据分析办法

李柏年吴礼斌主编张孔生丁华参编第1页第4章鉴别分析鉴别分析基本思想是根据已知类别样本所提供信息，总结出分类规律性，建立鉴别公式和鉴别准则，鉴别新样本点所属类型。本章介绍距离鉴别分析、Bayes鉴别分析极其MATLAB软件实现。4.1距离鉴别分析4.1.1鉴别分析概念在某些自然科学和社会科学研究中，研究对象用某种办法已划分为若干类型，当得到一种新样品数据（一般是多元），要确定该样品属于已知类型中哪一类，这样问题属于鉴别分析.第2页从统计数据分析角度，可概括为如下模型：设有k个总体,它们都是p元总体,其数量指标是1)若总体分布函数是已知，对于任一新样品数据,判断它来自哪一种总体。2)一般各个总体分布是未知，由从各个总体取得样本（训练样本）来估计。一般，先估计各个总体均值向量与协方差矩阵。第3页标准:1.从统计学角度，要求鉴别准则在某种准则下是最优，例如错判概率最小等。2.根据不一样鉴别准则，有不一样鉴别办法，这里主要介绍距离鉴别和Bayes鉴别4.1.2距离定义1.闵可夫斯基距离设有n维向量称绝对距离第4页称称为n维向量x，y之间闵可夫斯基距离，其中为常数。欧氏距离显然，当r=2和1时闵可夫斯基距离分别为欧氏距离和绝对距离.第5页(1)同一总体两个向量之间马氏距离其中

为总体协方差矩阵，一般取

为实对称正定矩阵.显然，当

为单位矩阵时马氏距离就是欧氏距离.设有n维向量,则称为n维向量x，y之间马氏距离.2.马氏距离马氏距离是由印度统计学家马哈拉诺比斯(PCMahalanobis)提出，由于马氏距离具有统计意义，在距离鉴别分析时经常应用马氏距离：(4.1.1)第6页(2)一种向量到一种总体马氏距离总体G均值向量为μ，协方差矩阵为Σ.则称为n维向量x与总体G马氏距离.MATLAB中有一种命令：d=mahal(Y,X)，计算X矩阵每一种点（行）至Y矩阵中每一种点（行）马氏距离。其中Y列数必须等于X列数，但它们行数能够不一样。X行数必须大于列数。输出d是距离向量。(4.1.2)第7页(3)两个总体之间马氏距离设有两个总体G1,G2，两个总体均值向量分别为，协方差矩阵相等，皆为Σ,则两个总体之间马氏距离为一般，在鉴别分析时不采取欧氏距离原因在于，该距离与量纲有关.例如平面上有A,B,C,D四个点，横坐标为代表重量（单位：kg），纵坐标代表长度（单位：cm），如下页图。

(4.1.3)第8页这时显然AB>CD假如目前长度用mm为单位，重量单位保持不变，于是A点坐标为(0,50)，B点坐标为(0,100)，此时计算线段长度为此时，AB<CD第9页这表白欧氏距离有一种缺陷，当向量分量是不一样量纲时欧氏距离大小竟然与指标单位有关.而马氏距离则与量纲无关.4.1.3两总体距离鉴别分析先考虑两个总体情况。设，为两个不一样p元已知总体，均值向量是，，协方差矩阵是，.设是一种待判样品，距离鉴别准则为(4.1.4)第10页即当到马氏距离不超出到马氏距离时，判来自；反之，判来自.由于马氏距离与总体协方差矩阵有关，因此利用马氏距离进行鉴别分析需要分别考虑两个总体协方差矩阵是否相等.1.两个总体协方差矩阵相等情况设有两个总体G1,G2，均值分别为,协方差矩阵相等为Σ。考虑样品x到两个总体马氏距离平方差：第11页其中，令于是距离鉴别准则为(4.1.6)第12页

由于总体均值、协方差矩阵一般是未知，数据资料来自两个总体训练样本，于是用样本均值、样本协方差矩阵替代总体均值与协方差.注意：若S1,S2分别为两个样本协方差矩阵，则在两个总体协方差矩阵相等时，总体协方差矩阵估计量其中n1,n2分别为两个样本容量.得到教材中鉴别法则：(4.1.11)(4.1.9)第13页matlab鉴别步骤：1.计算A、B两类均值向量与协方差阵;ma=mean(A),mb=mean(B),S1=cov(A),S2=cov(B)2.计算总体协方差矩阵其中n1,n2分别为两个样本容量.3.计算未知样本x到A,B两类马氏平方距离之差d=(x-ma)S-1(x-ma)’-(x-mb)S-1(x-mb)’4.若d<0,则x属于A类;若d>0,则x属于B类第14页上述公式能够化简为：W(x)=(ma-mb)S-1(x-(ma+mb)/2)’若W(x)>0，x属于G1;若W(x)<0，x属于G2注意：1.此处ma,mb都是行向量；2.当x是一种矩阵时，则用ones矩阵左乘(ma+mb)/2后来，方可与x相减.※Matlab中直接进行数据鉴别分析命令为classify，其调用格式class=classify(sample,training,group'type')第15页例4.1.1(1989年国际数学竞赛A题)蠓分类蠓是一种昆虫，分为很多类型，其中有一种名为Af,是能传输花粉益虫；另一种名为Apf,是会传输疾病害虫，这两种类型蠓在形态上十分相同，很难区分.现测得6只Apf和9只Af蠓虫触角长度和翅膀长度数据Apf：(1.14,1.78),(1.18,1.96),(1.20,1.86),(1.26,2.00),(1.28,2.00),(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).若两类蠓虫协方差矩阵相等，试鉴别下列三个蠓虫属于哪一类？(1.24,1.8)，(1.28,1.84)，（1.4，2.04）第16页解：假定两总体协方差相等，源程序如下：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:3W(i)=(x(i,:)-1/2*(m1+m2))*inv(s)*(m1-m2)';%计算鉴别函数值end;输出成果为：W=2.16401.35681.9802由鉴别准则(4.1.11)可知，三只蠓虫均属于Apf.第17页直接调用MATLAB鉴别分析命令classify。apf=[1.14,1.78;1.18,1.96;1.20,1.86;1.26,2.;1.28,2;1.30,1.96];%总体apfaf=[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];%总体aftraining=[apf;af];%合并两个总体形成训练集n1=size(apf,1);%总体apf中样本行数n2=size(af,1);%总体af中样本行数group=[ones(1,n1),2*ones(1,n2)];%apf中样本与af中样本类属x=[1.24,1.8;1.28,1.84;1.4,2.04];%输入原始待判数据即sampleclass=classify(x,training,group)%鉴别分析输出成果为：class=111由鉴别准则(4.1.11)可知，三只蠓虫均属于Apf.第18页2.两个总体协方差矩阵不相等样品到两个总体马氏距离平方分别为：令则鉴别准则：(4.1.13)第19页当两个总体协方差矩阵不等时，能够建立MATLAB鉴别法如下：例4.1.2对例4.1.1数据，假定两类总体协方差矩阵不相等，重新鉴别上述三个蠓虫类别.

解：程序如下：

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];%输入原始数据W=mahal(x,apf)-mahal(x,af)%计算鉴别函数第20页

输出成果为：

W=1.76113.88123.6468

由鉴别准则(4.1.17)可知，三个蠓虫均属于Af.3.两个总体协方差矩阵相等检查以上两个例题成果大相径庭，由此我们不禁要问究竟哪个成果更可靠？问题关键在于：两类蠓虫总体协方差矩阵是否相等？着手处理协方差矩阵检查.

检查统计量：第21页对给定，查卡方分布表得到临界值.若,则接收H0，不然回绝H0对于例4.1.1，应用检查程序如下：n1=6;n2=9;p=2;s=(5*s1+8*s2)/13;Q01=(n1-1)*(log(det(s))-log(det(s1))-p+trace(inv(s)*s1))；Q02=(n2-1)*(log(det(s))-log(det(s2))-p+trace(inv(s)*s2))；成果:Q01=2.5784，Q02=0.7418对，查自由度为3卡方分布chi2inv(0.05,3),得到临界值为：7.8147由于Q01<7.8147，Q02<7.8147，故以为两总体协方差矩阵相同。例4.1.1那种解法更合理.第22页4.1.4多种总体距离鉴别设有k个总体G1,G2,…,Gk，若鉴别某个体x属于哪个总体，则有如下办法：若存在某个正整数k0,使得mahal(y,Gk0)=min(mahal(y,Gi)),(i=1,2,…,k)则鉴别y属于第k0个总体.多种总体协方差矩阵是否相等检查(参照第二章第2.2.2节)第23页1.总体协方差矩阵相等时鉴别设有k个总体G1，G2，…，Gk，是取自总体Gj(j=1,2,…,k)训练样本，记于是未知样品到各总体鉴别函数为：其中鉴别准则为：若则x属于Gj0(4.1.21)第24页解：根据例2.2.3结论，能够以为三类总体协方差矩阵相等.A=[260 75 40 18 310 122 30 21 320 64 39 17;……260 135 39 29 280 40 37 17 250 117 36 16];G1=A(:,1:4);G2=A(:,5:8);G3=A(:,9:12);%三类总体数据x=[190673017;3151003519;240603718];%待判定数据m(1,:)=mean(G1);m(2,:)=mean(G2);m(3,:)=mean(G3);s1=cov(G1);s2=cov(G2);s3=cov(G3);s=19*(s1+s2+s3)/57;fori=1:3forj=1:3fork=1:3例4.1.3对例2.2.3表2.6中给出身体指标化验数据，对三个待判数(190,67,30,17),(315,100,35,19)，(240,60,37,18)进行鉴别归类。表2.6第25页

w(j,k)=(x(i,:)-1/2*(m(j,:)+m(k,:)))*inv(s)*(m(j,:)-m(k,:))';ifw(j,k)<0q=0;break;elseq=1;end;end;ifq==1y(i)=j;end;end;end;输出成果：y=132

由以上鉴别准则可知，三个待判数据（190,67,30,17），（315,100,35,19），（240,60,37,18）分别属于G1,G3和G2,和

第26页2.总体协方差矩阵不全相等时鉴别计算样品到各总体马氏距离平方：记鉴别准则为：若则鉴别x属于Gj0(4.1.27)第27页4.2鉴别准则评价

当一种鉴别准则提出后来，还要研究它优良性，即考查它误判概率.以训练样本为基础估计思想：若属于样品被误判为属于个数为个，属于样品被误判为属于个数为个，两类总体样品总数为，则误判概率估计为：

针对详细情况，一般采取回代法和交叉法进行误判概率估计。第28页（1）回代误判率设G1，G2为两个总体，X1,X2,…,Xm和Y1,Y2,…,Yn是分别来自G1，G2训练样本，以全体训练样本作为m+n个新样品，逐一代入已建立鉴别准则中鉴别其归属，这个过程称为回判。若属于G1样品被误判为属于G2个数为N1个，属于G2样品被误判为属于G1个数为N2个，则误判率估计为：(4.2.1)第29页（2）交叉误判率估计交叉误判率估计是每次剔除一种样品，利用其他m+n－1个训练样本建立鉴别准则再用所建立准则对删除样品进行鉴别。对训练样本中每个样品都做如上分析，以其误判百分比作为误判率，详细步骤如下：①从总体为G1训练样本开始，剔除其中一种样品，剩下m－1个样品与G2中所有样品建立鉴别函数；②用建立鉴别函数对剔除样品进行鉴别；第30页③反复步骤①，②，直到G1中所有样品依次被删除，又进行鉴别，其误判样品个数记为④对G2样品反复步骤①，②，③直到G2中所有样品依次被删除又进行鉴别，其误判样品个数记为于是交叉误判率估计为：(4.2.2)第31页例4.2.1

根据教材表4.1数据,鉴别两类总体协方差矩阵是否相等，然后用马氏距离鉴别未知地域类别，并计算回代误判率与交叉误判率.表4.1解：首先判断两组数据协方差是否相等.再建立鉴别准则，计算回代和交叉误判率，源程序如下：a=[503.1021.80332.30188.50…………769.9050.90605.0041.00];b=[89.709.50105.209.60…………1142.7030.80448.50334.20];x=[431.3047.20210.6014.40;1401.3047.20654.70350.701331.6057.00693.8020.40;279.9015.10118.505.10];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);第32页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))chi2inv(0.95,10)%验证两总体协方差矩阵相同fori=1:4D(i)=(x(i,:)-mean(a))*inv(s)*(x(i,:)-mean(a))'-(x(i,:)-mean(b))*inv(s)*(x(i,:)-mean(b))';end%由D成果可得：前三个属于第一类，最后一种属于第二类

fori=1:n1d11(i)=(a(i,:)-mean(a))*inv(s)*(a(i,:)-mean(a))'-(a(i,:)-mean(b))*inv(s)*(a(i,:)-mean(b))';endfori=1:n2d22(i)=(b(i,:)-mean(b))*inv(s)*(b(i,:)-mean(b))'-(b(i,:)-mean(a))*inv(s)*(b(i,:)-mean(a))';endn11=length(find(d11>0));n22=length(find(d22>0));p0=(n11+n22)/(n1+n2)%计算回代误判率第33页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(D11>0));N22=length(find(D22>0));p1=(N11+N22)/(n1+n2)%计算交叉误判率输出成果：p0=0.1923p1=0.2400

第34页4.3Bayes鉴别分析

贝叶斯公式是一种我们熟知公式

距离判别只要求知道总体数字特性，不涉及总体分布函数，当参数和协方差未知时，就用样本均值和协方差矩阵来估计。距离判别方法简单实用，但没有考虑到每个总体出现机会大小，即先验概率，没有考虑到错判损失。贝叶斯判别法正是为理解决这两个问题提出判别分析方法。第35页4.3.1两个总体Bayes鉴别1.一般讨论

考虑两个p元总体分别具有概率密度函数f1(x),f2(x)，设出现先验概率为：，且当取得新样品后,根据Bayes公式后验概率分别为(4.3.1)第36页因此，两个总体Bayes鉴别准则为2.两个正态总体Bayes鉴别（1）两个总体协方差矩阵相等情形设总体G1,G2协方差矩阵相等且为Σ，概率密度函数为：(4.3.2)第37页损失相等Bayes鉴别准则为其中基于两正态总体后验概率Bayes鉴别准则为其中第38页在实际问题中，有关先验概率，一般用下列两种方式选用:1）采取等概率选用，即2）按训练样本容量百分比选用，即第39页例4.3.1对例4.1.1数据，重新对上述三个蠓虫类别进行Bayes鉴别.（假设误判损失相等）解：第1步：能够验证两个总体服从二元正态分布；（第二章正态性检查，读者自证）第2步：检查两个总体协方差矩阵相等；第3步：估计两个总体先验概率,,这里按样本容量百分比选用.由于Apf与Af分别为6个与9个，故估计Apf类蠓虫先验概率，Af类蠓虫先验概率；第4步：利用MATLAB软件计算：第40页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:3

w1(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类']);

else

disp(['第',num2str(i),'个蠓虫属于Af类']);

end;end;输出成果：

第1个蠓虫属于Apf类第2个蠓虫属于Apf类第3个蠓虫属于Apf类第41页（2）两个总体协方差矩阵不相等情形设总体协方差矩阵不相等分别为Σ1,Σ2概率密度函数为：则基于两正态总体误判损失相等Bayes鉴别准则其中第42页

例4.3.2对破产公司搜集它们在破产前两年年度财务数据，对财务良好公司也搜集同一时间数据.数据包括四个变量:钞票流量/总债务，净收益/总资产，流动资产/流动债务，以及流动资产/净销售额，数据如表4.2所示.假定两总体G1,G2均服从四元正态分布，在误判损失相等且先验概率按百分比分派条件下，看待判样本进行bayes鉴别.表4.2解：第1步：检查两个总体协方差矩阵相等；源程序如下：第43页A=[-0.45-0.411.090.45 0.51 0.10 2.49 0.54……-0.13-0.14 1.42 0.44 0.17 0.07 1.80 0.52]x=[-0.23-0.30 0.330.18;0.15 0.05 2.17 0.55-0.28-0.23 1.190.66;0.48 0.09 1.24 0.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));ifQ1<chi2inv(0.95,p*(p+1)/2)&&Q2<chi2inv(0.95,p*(p+1)/2)disp('两组数据协方差相等');elsedisp('两组数据协方差不全相等');end;输出成果：两组数据协方差不全相等第44页第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个属于非破产公司第45页4.3.2多种总体Bayes鉴别设有k个总体G1,G2,…,Gk概率密度为fj(x)各总体出现先验概率为1.一般讨论当出现样品时,总体后验概率第46页Bayes鉴别准则为:若则判样本注：当达成最大后验概率不止一种时，可判为达成最大后验概率总体任何一种.2.多种正态总体Bayes鉴别（1）当时，设

第47页线性鉴别函数为其中基于误判损失相等Bayes鉴别准则为基于后验概率Bayes鉴别准则为其中第48页在实际问题中，由于未知，各总体训练样本均值（2）当不全相等时,设则基于后验概率Bayes鉴别准则为其中未知，估计.第49页例4.3.3.某医院利用心电图检测来对人群进行划分，数据见表.“g=1”表达健康人，“g=2”表达积极脉硬化患者，“g=3”表达冠心病患者，X1

,X2表达测得心电图中表白心脏功能两项不有关指标.某受试者心电图该两项指标数据分别为380.20，9.08.设先验概率按百分比分派，进行bayes鉴别，判定其归属.表4.324人心电图数据编号X1X2编号X1X2123456789101112261.01185.39249.58137.13231.34231.38260.25259.51273.84303.59231.03308.907.365.996.114.358.798.5310.029.798.798.536.158.49111111111112131415161718192021222324258.69355.54476.69316.12274.57409.42330.34331.47352.50347.31189.59380.207.169.4311.328.179.6710.499.6113.7211.0011.195.469.0822222233333待判第50页解：A=[261.01 7.36185.395.99……189.59 5.46]x=[380.20 9.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*p^2+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%计算统计量观测值第51页C=chi2inv(0.95,f)ifT<chi2inv(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组

第52页4.3.3平均误判率Byaes鉴别有效性能够通过平均误判率来确定。这里仅对两个正态总体,且协方差矩阵相等情况下研究平均误判率计算.设总体，其先验概率，两个总体马氏平方距离记为第53页则基于误判损失相等时平均误判率为其中为标准正态分布函数.从（4.3.11）式知，当总体马氏平方距离越大，即两总体分离程度越大时，平均误判概率最小.推广到一般情况也成立.（4.3.11）第54页例4.3.42023年全国部分地域城镇居民人均年家收入情况见表.按四种指标分为二类，用bayes鉴别判定青海、广东两省区属于哪一类，并用回代法和交叉法对误判率进行估计(假定误判损失相等).表4.4解：第1步，检查三个总体协方差矩阵相等；A=[18738.96 778.36 452.75 7707.87……9422.22 938.15 141.75 1976.49];x=[8595.48763.07 50.17 3458.6315188.392405.92701.25 3382.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样本数第55页n=n1+n2+n3;%三个总体合并样本数k=3;p=4;f=p*(p+1)*(k-1)/2;d=(2*p^2+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)ifT<chi2inv(0.95,f)disp('三组数据协方差相等');elsedisp('三组数据协方差不全相等');end输出成果：三组数据协方差相等第56页第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;第57页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,:))&j~=1

n11=n11+1;endendendfori=1:n2

w2(i,1)=m1*inv(s)*G2(i,:)'-1/2*m1*inv(s)*m1'+log(p1);w2(i,2)=m2*inv(s)*G2(i,:)'-1/2*m2*inv(s)*m2'+log(p2);w2(i,3)=m3*inv(s)*G2(i,:)'-1/2*m3*inv(s)*m3'+log(p3);%计算鉴别函数

forj=1:3ifw2(i,j)==max(w2(i,:))&j~=2

n22=n22+1;endendend

第58页fori=1:n3

w3(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,:))&j~=3

n33=n33+1;endendendp00=(n11+n22+n33)/(n1+n2+n3)输出成果：p00=0第4步，计算交叉误判率.N11=0;N22=0;N33=0;

第59页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*in