MATLAB數(shù)據(jù)分析方法-判別分析

MATLAB數(shù)據(jù)分析方法

第4章判別分析判別分析的基本思想是根據(jù)已知類別的樣本所提供的信息，總結出分類的規(guī)律性，建立判別公式和判別準則，判別新的樣本點所屬類型。本章介紹距離判別分析、Bayes判別分析極其MATLAB軟件的實現(xiàn)。

4.1距離判別分析4.1.1判別分析的概念在一些自然科學和社會科學的研究中，研究對象用某種方法已劃分為若干類型，當?shù)玫降囊粋€新樣品數(shù)據(jù)（通常是多元的），要確定該樣品屬于已知類型中的哪一類，這樣的問題屬于判別分析.從統(tǒng)計數(shù)據(jù)分析的角度，可概括為如下模型：設有k個總體,它們都是p元總體,其數(shù)量指標是

1)若總體的分布函數(shù)是已知，對于任一新樣品數(shù)據(jù),判斷它來自哪一個總體。

2)通常各個總體的分布是未知的，由從各個總體取得的樣本（訓練樣本）來估計。一般，先估計各個總體的均值向量與協(xié)方差矩陣。原則:1.從統(tǒng)計學的角度，要求判別準則在某種準則下是最優(yōu)的，例如錯判的概率最小等。

2.根據(jù)不同的判別準則，有不同的判別方法，這里主要介紹距離判別和Bayes判別

4.1.2距離的定義

1.閔可夫斯基距離設有n維向量稱絕對距離稱稱為n維向量x，y之間的閔可夫斯基距離，其中為常數(shù)。歐氏距離顯然，當r=2和1時閔可夫斯基距離分別為歐氏距離和絕對距離.(1)同一總體的兩個向量之間的馬氏距離其中

為總體協(xié)方差矩陣，通常取

為實對稱正定矩陣.顯然，當

為單位矩陣時馬氏距離就是歐氏距離.設有n維向量,則稱為n維向量x，y之間的馬氏距離.2.馬氏距離馬氏距離是由印度統(tǒng)計學家馬哈拉諾比斯(PCMahalanobis)提出的，由于馬氏距離具有統(tǒng)計意義，在距離判別分析時經(jīng)常應用馬氏距離：(4.1.1)(2)一個向量到一個總體的馬氏距離總體G

的均值向量為μ，協(xié)方差矩陣為Σ.則稱為n維向量x與總體G的馬氏距離.

MATLAB中有一個命令：d=mahal(Y,X)，計算X矩陣每一個點（行）至Y矩陣中每一個點（行）的馬氏距離。其中Y的列數(shù)必須等于X的列數(shù)，但它們的行數(shù)可以不同。X的行數(shù)必須大于列數(shù)。輸出d是距離向量。

(4.1.2)(3)兩個總體之間的馬氏距離設有兩個總體G1,G2，兩個總體的均值向量分別為

，協(xié)方差矩陣相等，皆為Σ,則兩個總體之間的馬氏距離為通常，在判別分析時不采用歐氏距離的原因在于，該距離與量綱有關.例如平面上有A,B,C,D四個點，橫坐標為代表重量（單位：kg），縱坐標代表長度（單位：cm），如下頁圖。

(4.1.3)這時顯然AB>CD如果現(xiàn)在長度用mm為單位，重量的單位保持不變，于是A點的坐標為(0,50)，B點的坐標為(0,100)，此時計算線段的長度為此時，AB<CD這表明歐氏距離有一個缺陷，當向量的分量是不同的量綱時歐氏距離的大小竟然與指標的單位有關.而馬氏距離則與量綱無關.

4.1.3兩總體的距離判別分析先考慮兩個總體的情況。設，為兩個不同的p元已知總體，的均值向量是，，的協(xié)方差矩陣是，.設是一個待判樣品，距離判別準則為(4.1.4)即當?shù)降鸟R氏距離不超過到的馬氏距離時，判來自；反之，判來自.由于馬氏距離與總體的協(xié)方差矩陣有關，所以利用馬氏距離進行判別分析需要分別考慮兩個總體的協(xié)方差矩陣是否相等.1.兩個總體協(xié)方差矩陣相等的情況設有兩個總體G1,G2，均值分別為,協(xié)方差矩陣相等為Σ?？紤]樣品x到兩個總體的馬氏距離平方差：其中，令于是距離判別準則為(4.1.6)

由于總體的均值、協(xié)方差矩陣通常是未知的，數(shù)據(jù)資料來自兩個總體的訓練樣本，于是用樣本的均值、樣本的協(xié)方差矩陣代替總體的均值與協(xié)方差.注意：若S1,S2分別為兩個樣本的協(xié)方差矩陣，則在兩個總體協(xié)方差矩陣相等時，總體的協(xié)方差矩陣估計量其中n1,n2分別為兩個樣本的容量.得到教材中判別法則：(4.1.11)

(4.1.9)matlab判別步驟：

1.計算A、B兩類的均值向量與協(xié)方差陣;ma=mean(A),mb=mean(B),S1=cov(A),S2=cov(B)2.計算總體的協(xié)方差矩陣其中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類上述公式可以化簡為：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中直接進行數(shù)據(jù)的判別分析命令為classify，其調用格式class=classify(sample,training,group'type')例4.1.1(1989年國際數(shù)學競賽A題)蠓的分類蠓是一種昆蟲，分為很多類型，其中有一種名為Af,是能傳播花粉的益蟲；另一種名為Apf,是會傳播疾病的害蟲，這兩種類型的蠓在形態(tài)上十分相似，很難區(qū)別.現(xiàn)測得6只Apf和9只Af蠓蟲的觸角長度和翅膀長度數(shù)據(jù)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).

若兩類蠓蟲協(xié)方差矩陣相等，試判別以下的三個蠓蟲屬于哪一類？(1.24,1.8)，(1.28,1.84)，（1.4，2.04）解：假定兩總體的協(xié)方差相等，源程序如下：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];%輸入原始數(shù)據(jù)m1=mean(apf);

m2=mean(af);s1=cov(apf);s2=cov(af);s=(5*s1+8*s2)/13;%計算樣本均值與協(xié)方差矩陣fori=1:3W(i)=(x(i,:)-1/2*(m1+m2))*inv(s)*(m1-m2)';%計算判別函數(shù)值

end;輸出結果為：W=2.16401.35681.9802由判別準則(4.1.11)可知，三只蠓蟲均屬于Apf.直接調用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中樣本的行數(shù)n2=size(af,1);

%總體af中樣本的行數(shù)group=[ones(1,n1),2*ones(1,n2)];

%apf中樣本與af中樣本類屬x=[1.24,1.8;1.28,1.84;1.4,2.04];

%輸入原始待判數(shù)據(jù)即sampleclass=classify(x,training,group)

%判別分析輸出結果為：class=111由判別準則(4.1.11)可知，三只蠓蟲均屬于Apf.2.兩個總體協(xié)方差矩陣不相等樣品到兩個總體的馬氏距離平方分別為：令則判別準則：(4.1.13)當兩個總體的協(xié)方差矩陣不等時，可以建立MATLAB的判別法如下：例4.1.2對例4.1.1的數(shù)據(jù)，假定兩類總體的協(xié)方差矩陣不相等，重新判別上述三個蠓蟲的類別.

解：程序如下：

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];%輸入原始數(shù)據(jù)

W=mahal(x,apf)-mahal(x,af)%計算判別函數(shù)

輸出結果為：

W=1.76113.88123.6468

由判別準則(4.1.17)可知，三個蠓蟲均屬于Af.3.兩個總體協(xié)方差矩陣相等的檢驗

以上兩個例題的結果大相徑庭，由此我們不禁要問究竟哪個結果更可靠？問題的關鍵在于：兩類蠓蟲總體的協(xié)方差矩陣是否相等？著手解決協(xié)方差矩陣的檢驗.

檢驗統(tǒng)計量：對給定的，查卡方分布表得到臨界值.若,則接受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，故認為兩總體協(xié)方差矩陣相同。例4.1.1的那種解法更合理.

4.1.4多個總體的距離判別設有k個總體G1,G2,…,Gk，若判別某個體x屬于哪個總體，則有如下方法：若存在某個正整數(shù)k0,使得mahal(y,Gk0)=min(mahal(y,Gi)),(i=1,2,…,k)則判別y屬于第k0個總體.多個總體協(xié)方差矩陣是否相等的檢驗(參考第二章第2.2.2節(jié))1.總體協(xié)方差矩陣相等時的判別設有k個總體G1，G2，…，Gk，是取自總體Gj

(j=1,2,…,k)的訓練樣本，記于是未知樣品到各總體的判別函數(shù)為：其中判別準則為：若則x屬于Gj0(4.1.21)解：根據(jù)例2.2.3的結論，可以認為三類總體協(xié)方差矩陣相等.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);%三類總體數(shù)據(jù)x=[190673017;3151003519;240603718];%待判定的數(shù)據(jù)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:3

fork=1:3例4.1.3對例2.2.3表2.6中給出的身體指標化驗數(shù)據(jù)，對三個待判數(shù)(190,67,30,17),(315,100,35,19)，(240,60,37,18)進行判別歸類。

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

由以上判別準則可知，三個待判數(shù)據(jù)（190,67,30,17），（315,100,35,19），（240,60,37,18）分別屬于G1,G3和G2,和

2.總體協(xié)方差矩陣不全相等時的判別計算樣品到各總體馬氏距離平方：記判別準則為：若則判別x屬于Gj0(4.1.27)4.2判別準則的評價

當一個判別準則提出以后，還要研究它的優(yōu)良性，即考察它的誤判概率.以訓練樣本為基礎的估計思想：若屬于的樣品被誤判為屬于的個數(shù)為個，屬于的樣品被誤判為屬于的個數(shù)為個，兩類總體的樣品總數(shù)為，則誤判概率的估計為：

針對具體情況，通常采用回代法和交叉法進行誤判概率的估計。（1）回代誤判率設G1，G2為兩個總體，X1,X2,…,Xm和Y1,Y2,…,Yn是分別來自G1，G2的訓練樣本，以全體訓練樣本作為m+n個新樣品，逐個代入已建立的判別準則中判別其歸屬，這個過程稱為回判。若屬于G1的樣品被誤判為屬于G2的個數(shù)為N1個，屬于G2的樣品被誤判為屬于G1的個數(shù)為N2個，則誤判率估計為：

(4.2.1)（2）交叉誤判率估計交叉誤判率估計是每次剔除一個樣品，利用其余的m+n－1個訓練樣本建立判別準則再用所建立的準則對刪除的樣品進行判別。對訓練樣本中每個樣品都做如上分析，以其誤判的比例作為誤判率，具體步驟如下：①

從總體為G1的訓練樣本開始，剔除其中一個樣品，剩余的m－1個樣品與G2中的全部樣品建立判別函數(shù)；②

用建立的判別函數(shù)對剔除的樣品進行判別；③重復步驟①，②，直到G1中的全部樣品依次被刪除，又進行判別，其誤判的樣品個數(shù)記為④對G2的樣品重復步驟①，②，③直到G2中的全部樣品依次被刪除又進行判別，其誤判的樣品個數(shù)記為于是交叉誤判率估計為：

(4.2.2)例4.2.1

根據(jù)教材表4.1數(shù)據(jù),判別兩類總體的協(xié)方差矩陣是否相等，然后用馬氏距離判別未知地區(qū)的類別，并計算回代誤判率與交叉誤判率.解：首先判斷兩組數(shù)據(jù)協(xié)方差是否相等.再建立判別準則，計算回代和交叉誤判率，源程序如下：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);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)

%驗證兩總體的協(xié)方差矩陣相同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)%計算回代誤判率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

4.3Bayes判別分析

貝葉斯公式是一個我們熟知的公式

距離判別只要求知道總體的數(shù)字特征，不涉及總體的分布函數(shù)，當參數(shù)和協(xié)方差未知時，就用樣本的均值和協(xié)方差矩陣來估計。距離判別方法簡單實用，但沒有考慮到每個總體出現(xiàn)的機會大小，即先驗概率，沒有考慮到錯判的損失。貝葉斯判別法正是為了解決這兩個問題提出的判別分析方法。

4.3.1兩個總體的Bayes判別1.一般討論

考慮兩個p元總體分別具有概率密度函數(shù)f1(x),f2(x)，設出現(xiàn)的先驗概率為：，且當取得新樣品后,根據(jù)Bayes公式的后驗概率分別為

(4.3.1)因此，兩個總體的Bayes判別準則為2.兩個正態(tài)總體的Bayes判別（1）兩個總體協(xié)方差矩陣相等的情形設總體G1,G2的協(xié)方差矩陣相等且為Σ，概率密度函數(shù)為：(4.3.2)損失相等的Bayes判別準則為其中基于兩正態(tài)總體后驗概率的Bayes判別準則為其中在實際問題中，關于先驗概率，通常用下列兩種方式選取:1）采用等概率選取，即2）按訓練樣本的容量的比例選取，即例4.3.1

對例4.1.1的數(shù)據(jù)，重新對上述三個蠓蟲的類別進行Bayes判別.（假設誤判損失相等）第1步：可以驗證兩個總體服從二元正態(tài)分第2步：檢驗兩個總體的協(xié)方差矩陣相等；第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: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類（2）兩個總體協(xié)方差矩陣不相等的情形設總體的協(xié)方差矩陣不相等分別為Σ1,Σ2概率密度函數(shù)為：則基于兩正態(tài)總體誤判損失相等的Bayes判別準則其中

例4.3.2對破產的企業(yè)收集它們在破產前兩年的年度財務數(shù)據(jù)，對財務良好的企業(yè)也收集同一時間的數(shù)據(jù).數(shù)據(jù)涉及四個變量:現(xiàn)金流量/總債務，

凈收益/總資產，

流動資產/流動債務，以及

流動資產/凈銷售額，數(shù)據(jù)如表4.2

所示.假定兩總體G1,G2均服從四元正態(tài)分布，在誤判損失相等且先驗概率按比例分配的條件下，對待判樣本進行bayes判別.

解：第1步：檢驗兩個總體的協(xié)方差矩陣相等；源程序如下：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);%二類總體數(shù)據(jù)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('兩組數(shù)據(jù)協(xié)方差相等');elsedisp('兩組數(shù)據(jù)協(xié)方差不全相等');end;輸出結果：

兩組數(shù)據(jù)協(xié)方差不全相等第2步：根據(jù)第1步結論，構造判別函數(shù)，得出判結果.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),'個屬于破產企業(yè)']);elsedisp(['第',num2str(i),'個屬于非破產企業(yè)']);end;end;輸出結果：第1個屬于破產企業(yè)第2個屬于非破產企業(yè)第3個屬于破產企業(yè)第4個屬于非破產企業(yè)4.3.2多個總體的Bayes判別設有k個總體G1,G2,…,Gk的概率密度為fj(x)各總體出現(xiàn)的先驗概率為

1.一般討論當出現(xiàn)樣品時,總體的后驗概率Bayes判別準則為:若則判樣本注：當達到最大后驗概率的不止一個時，可判為達到最大后驗概率的總體的任何一個.2.多個正態(tài)總體的Bayes判別（1）當時，設

線性判別函數(shù)為其中基于誤判損失相等的Bayes判別準則為基于后驗概率的Bayes判別準則為其中在實際問題中，由于未知，各總體的訓練樣本均值（2）當不全相等時,設則基于后驗概率的Bayes判別準則為其中未知，估計.例4.3.3.某醫(yī)院利用心電圖檢測來對人群進行劃分，數(shù)據(jù)見表.“g=1”表示健康人，“g=2”表示主動脈硬化患者，“g=3”表示冠心病患者，X1

,X2表示測得的心電圖中表明心臟功能的兩項不相關的指標.某受試者心電圖該兩項指標的數(shù)據(jù)分別為380.20，9.08.設先驗概率按比例分配，進行bayes判別，判定其歸屬.表4.324人心電圖數(shù)據(jù)編號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待判解：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,:);%三類總體數(shù)據(jù)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);%計算協(xié)方差陣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%計算統(tǒng)計量觀測值C=chi2inv(0.95,f)ifT<chi2inv(0.95,f)

disp('三組數(shù)據(jù)協(xié)方差相等');else

disp('三組數(shù)據(jù)協(xié)方差不全相等');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;輸出結果:三組數(shù)據(jù)協(xié)方差相等屬于第2組

4.3.3平均誤判率Byaes判別的有效性可以通過平均誤判率來確定。這里僅對兩個正態(tài)總體,且協(xié)方差矩陣相等的情況下研究平均誤判率的計算.設總體，其先驗概率，兩個總體的馬氏平方距離記為則基于誤判損失相等時的平均誤判率為其中為標準正態(tài)分布函數(shù).從（4.3.11）式知，當總體的馬氏平方距離越大，即兩總體的分離程度越大時，平均誤判概率最小.推廣到一般情況也成立.（4.3.11）例4.3.42008年全國部分地區(qū)城鎮(zhèn)居民人均年家收入情況見表.按四種指標分為二類，用bayes判別判定青海、廣東兩省區(qū)屬于哪一類，并用回代法和交叉法對誤判率進行估計(假定誤判損失相等).解：第1步，檢驗三個總體的協(xié)方差矩陣相等；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.39

2405.92

701.25 3382.95];

%待判樣品G1=A(1:2,:);G2=A(3:8,:);G3=A(9:27,:);%輸入三類總體數(shù)據(jù)n1=size(G1,1);%總體G1的樣本數(shù)n2=size(G2,1);%總體G2的樣本數(shù)n3=size(G3,1);%總體G3的樣本數(shù)n=n1+n2+n3;

%三個總體合并的樣本數(shù)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);%計算協(xié)方差陣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%計算統(tǒng)計量觀測值C=chi2inv(0.95,f)ifT<chi2inv(0.95,f)disp('三組數(shù)據(jù)協(xié)方差相等');elsedisp('三組數(shù)據(jù)協(xié)方差不全相等');end輸出結果：三組數(shù)據(jù)協(xié)方差相等第2步，根據(jù)第1步結論，構造判別函數(shù)，得出判別結果.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);%計算判別函數(shù)

forj=1:3

ifw(j)==max(w)

disp(['待判樣品屬于第',num2str(j),'類城市']);

end

endend輸出結果：待判樣品屬于第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);%計算判別函數(shù)

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);%計算判別函數(shù)

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

n22=n22+1;endendend

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);%計算判別函數(shù)

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;

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