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

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

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

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

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

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

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

4.1.2距離的定義

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

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

為實(shí)對(duì)稱正定矩陣.顯然，當(dāng)

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

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

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

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

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

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

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

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

(4.1.9)matlab判別步驟：

1.計(jì)算A、B兩類的均值向量與協(xié)方差陣;ma=mean(A),mb=mean(B),S1=cov(A),S2=cov(B)2.計(jì)算總體的協(xié)方差矩陣其中n1,n2分別為兩個(gè)樣本的容量.3.計(jì)算未知樣本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類上述公式可以化簡(jiǎn)為：W(x)=(ma-mb)S-1(x-(ma+mb)/2)’若W(x)>0，x屬于G1;若W(x)<0，x屬于G2注意：1.此處ma,mb都是行向量；2.當(dāng)x是一個(gè)矩陣時(shí)，則用ones矩陣左乘(ma+mb)/2以后，方可與x相減.※Matlab中直接進(jìn)行數(shù)據(jù)的判別分析命令為classify，其調(diào)用格式class=classify(sample,training,group'type')例4.1.1(1989年國(guó)際數(shù)學(xué)競(jìng)賽A題)蠓的分類蠓是一種昆蟲(chóng)，分為很多類型，其中有一種名為Af,是能傳播花粉的益蟲(chóng)；另一種名為Apf,是會(huì)傳播疾病的害蟲(chóng)，這兩種類型的蠓在形態(tài)上十分相似，很難區(qū)別.現(xiàn)測(cè)得6只Apf和9只Af蠓蟲(chóng)的觸角長(zhǎng)度和翅膀長(zhǎng)度數(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).

若兩類蠓蟲(chóng)協(xié)方差矩陣相等，試判別以下的三個(gè)蠓蟲(chóng)屬于哪一類？(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;%計(jì)算樣本均值與協(xié)方差矩陣fori=1:3W(i)=(x(i,:)-1/2*(m1+m2))*inv(s)*(m1-m2)';%計(jì)算判別函數(shù)值

end;輸出結(jié)果為：W=2.16401.35681.9802由判別準(zhǔn)則(4.1.11)可知，三只蠓蟲(chóng)均屬于Apf.直接調(diào)用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];

%合并兩個(gè)總體形成訓(xùn)練集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)

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

解：程序如下：

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

輸出結(jié)果為：

W=1.76113.88123.6468

由判別準(zhǔn)則(4.1.17)可知，三個(gè)蠓蟲(chóng)均屬于Af.3.兩個(gè)總體協(xié)方差矩陣相等的檢驗(yàn)

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

檢驗(yàn)統(tǒng)計(jì)量：對(duì)給定的，查卡方分布表得到臨界值.若,則接受H0，否則拒絕H0對(duì)于例4.1.1

，應(yīng)用檢驗(yàn)程序如下：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))；結(jié)果:Q01=2.5784，Q02=0.7418對(duì)，查自由度為3的卡方分布chi2inv(0.05,3),得到臨界值為：7.8147由于Q01<7.8147，Q02<7.8147，故認(rèn)為兩總體協(xié)方差矩陣相同。例4.1.1的那種解法更合理.

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

(j=1,2,…,k)的訓(xùn)練樣本，記于是未知樣品到各總體的判別函數(shù)為：其中判別準(zhǔn)則為：若則x屬于Gj0(4.1.21)解：根據(jù)例2.2.3的結(jié)論，可以認(rèn)為三類總體協(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對(duì)例2.2.3表2.6中給出的身體指標(biāo)化驗(yàn)數(shù)據(jù)，對(duì)三個(gè)待判數(shù)(190,67,30,17),(315,100,35,19)，(240,60,37,18)進(jìn)行判別歸類。

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;輸出結(jié)果：y=132

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

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

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

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

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

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

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

(4.2.2)例4.2.1

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

%驗(yàn)證兩總體的協(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結(jié)果可得：前三個(gè)屬于第一類，最后一個(gè)屬于第二類

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)%計(jì)算回代誤判率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)%計(jì)算交叉誤判率輸出結(jié)果：p0=0.1923p1=0.2400

4.3Bayes判別分析

貝葉斯公式是一個(gè)我們熟知的公式

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

4.3.1兩個(gè)總體的Bayes判別1.一般討論

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

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

對(duì)例4.1.1的數(shù)據(jù)，重新對(duì)上述三個(gè)蠓蟲(chóng)的類別進(jìn)行Bayes判別.（假設(shè)誤判損失相等）第1步：可以驗(yàn)證兩個(gè)總體服從二元正態(tài)分第2步：檢驗(yàn)兩個(gè)總體的協(xié)方差矩陣相等；第3步：估計(jì)兩個(gè)總體的先驗(yàn)概率,,這里按樣本容量的比例選取.由于Apf與Af分別為6個(gè)與9個(gè)，故估計(jì)Apf類蠓蟲(chóng)的先驗(yàn)概率，Af類蠓蟲(chóng)的先驗(yàn)概率

；第4步：利用MATLAB軟件計(jì)算：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),'個(gè)蠓蟲(chóng)屬于Apf類']);

else

disp(['第',num2str(i),'個(gè)蠓蟲(chóng)屬于Af類']);

end;end;

輸出結(jié)果：

第1個(gè)蠓蟲(chóng)屬于Apf類

第2個(gè)蠓蟲(chóng)屬于Apf類

第3個(gè)蠓蟲(chóng)屬于Apf類（2）兩個(gè)總體協(xié)方差矩陣不相等的情形設(shè)總體的協(xié)方差矩陣不相等分別為Σ1,Σ2概率密度函數(shù)為：則基于兩正態(tài)總體誤判損失相等的Bayes判別準(zhǔn)則其中

例4.3.2對(duì)破產(chǎn)的企業(yè)收集它們?cè)谄飘a(chǎn)前兩年的年度財(cái)務(wù)數(shù)據(jù)，對(duì)財(cái)務(wù)良好的企業(yè)也收集同一時(shí)間的數(shù)據(jù).數(shù)據(jù)涉及四個(gè)變量:現(xiàn)金流量/總債務(wù)，

凈收益/總資產(chǎn)，

流動(dòng)資產(chǎn)/流動(dòng)債務(wù)，以及

流動(dòng)資產(chǎn)/凈銷售額，數(shù)據(jù)如表4.2

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

解：第1步：檢驗(yàn)兩個(gè)總體的協(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;輸出結(jié)果：

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

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

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

,X2表示測(cè)得的心電圖中表明心臟功能的兩項(xiàng)不相關(guān)的指標(biāo).某受試者心電圖該兩項(xiàng)指標(biāo)的數(shù)據(jù)分別為380.20，9.08.設(shè)先驗(yàn)概率按比例分配，進(jìn)行bayes判別，判定其歸屬.表4.324人心電圖數(shù)據(jù)編號(hào)X1X2編號(hào)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);%計(jì)算協(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%計(jì)算統(tǒng)計(jì)量觀測(cè)值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;輸出結(jié)果:三組數(shù)據(jù)協(xié)方差相等屬于第2組

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

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

forj=1:3

ifw(j)==max(w)

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

end

endend輸出結(jié)果：待判樣品屬于第3類城市

待判樣品屬于第2類城市第3步，計(jì)算回代誤判率.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);%計(jì)算判別函數(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);%計(jì)算判別函數(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);%計(jì)算判別函數(shù)

forj=1:3ifw3(i,j)==max(w3(i,:))&j~=3

n33=n33+1;endendendp00=(n11+n22+n33)/(n1+n2+n3)輸出結(jié)果：p00=

0第4步，計(jì)算交叉誤判率.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