主編：費宇中國人民大學(xué)出版社第8章判別分析2024/7/31

2判別分析是在已知樣品所有可能分類的前提下，將給定的新樣品按照某種分類準則判入其中某個類中的一種多元統(tǒng)計方法．例如：根據(jù)患者的各項檢查指標來判斷該病人屬于哪類病癥；根據(jù)某地氣象的記錄資料來判別(預(yù)報)未來幾天的天氣狀況；根據(jù)某地相關(guān)經(jīng)濟指標判斷該地區(qū)屬于哪一種經(jīng)濟類型地區(qū)．2024/7/318.1距離判別距離判別簡介

、

兩個及多個總體的距離判別8.2

Fisher判別兩個及多個總體的Fisher判別8.3

Bayes判別兩個及多個總體的Bayes判別8.4二次判別8.5案例分析與R實現(xiàn)本章主要內(nèi)容：38.1距離判別思想：根據(jù)一個樣品與各個類別距離的遠近對

其所屬類別進行判定．這里的距離通常使用馬氏距離：4

2024/7/315

2024/7/316

2024/7/317

2024/7/31例8.1

(冠心病例指標判斷)2024/7/31

8

測定了50-59歲冠心病人15例和正常人15例的舒張壓x1和膽固醇指標x2.試據(jù)此數(shù)據(jù)做距離判別分析．今若測得兩個患者的兩項指標分別為(90,160)和(85,155)，對他們?nèi)绾芜M行判斷？冠心病人正常人組別x1x2組別x1x2174200294172110014421001181110150270152170274280172196212280190180158270142180172280107110014028012411002302801941100220278152190239270190111015528010411001552809419614028413211002302701402024/7/31

9>setwd("C:/data")#設(shè)定工作路徑>d8.1<-read.csv("exam8.1.csv",header=T)#將exam8.1數(shù)據(jù)讀入>A1=d8.1[1:15,3:4];A2=d8.1[16:30,3:4]#將兩個總體樣本分開>mu1=apply(A1,2,mean);mu2=apply(A2,2,mean);mu1;mu2#計算兩個總體樣本的均值x1x293.73333187.93333x1x279.73333145.53333>S1=var(A1);S2=var(A2);S1;S2#計算兩個總體樣的協(xié)方差矩陣x1x2x1151.3524-221.5905x2-221.59051899.3524x1x2x172.49524-47.5619x2-47.561901078.69522024/7/31

10>W2unequal=function(x,mu1,mu2,S1,S2){mahalanobis(x,mu2,S2)-mahalanobis(x,mu1,S1)}>x1=c(90,160);x2=c(85,155)>W2unequal(x1,mu1,mu2,S1,S2)[1]1.082728#將x1代入值為正，判斷該點屬于G1（冠心病人）>W2unequal(x2,mu1,mu2,S1,S2)[1]-1.289129#將x2代入，值為負，判斷該點屬于G2（正常人）作圖進行直觀理解：>plot(d8.1[,3:4],type="n")>points(d8.1[1:15,3:4],pch=16)#用實心小圓點標出冠心病樣品點>points(d8.1[16:30,3:4],pch=21)#用空心小圓點標出正常人樣品點>points(90,160,pch=17)#用實心三角點標出第1位患者樣品點>points(85,155,pch=24)#用空心三角點標出第2位患者樣品點2024/7/31

11圖8－2

兩個總體樣品點以及兩位被判患者的樣品點多個總體的距離判別

2024/7/31

128.2

Fisher判別Fisher于1936年在植物分類研究中提出了該判別法．主要思想是通過將多維數(shù)據(jù)投影到一維直線上，使得同一類別(總體)中的數(shù)據(jù)在該直線上盡量靠攏，不同類別(總體)的數(shù)據(jù)盡可能分開，然后再利用前面的距離判別法來建立判別準則．這種投影直線（判別函數(shù)）可能有一條或多條．常用方法：線性判別法；非線性判別法；典型判別法等．這里主要介紹線性判別法．2024/7/31

135.2.1兩總體的Fisher判別2024/7/31

14

兩總體Fisher判別的思想是將二維空間中的點投影到一條適當?shù)闹本€

y

上，使得兩個總體G1和

G2中的點在直線y上的投影點盡可能分開，而同一總體在y上的投影點盡可能靠攏，再利用上面的距離判別法來建立判別準則．

示意圖：總體G1（大圓點）和

G2

（小圓點）中的點在直線y上的投影點滿足上面要求．再根據(jù)新樣品點x在

y上的投影點

y=aTx靠近μ1y（在μy右側(cè)）還是靠近μ2y（在μy左側(cè)）來判斷新樣品點x的歸屬．2024/7/31

15圖8-3兩總體Fisher判別示意圖兩總體的Fisher判別準則：2024/7/31

16

8.2.2

多總體Fisher判別多總體情形，通常要選取多條投影直線，即選取多個判別函數(shù)

來進行判別．設(shè)有k個總體G1,G2,···,

Gk，它們有共同的協(xié)方差陣Σ，均值分別為μ1,μ2

,···,μk.令問題為選擇a,使的投影變差比值達到最大:2024/7/31

178.2.2

多總體Fisher判別定理8.2設(shè)為的s個非零特征根，為相應(yīng)的特征向量且滿足

，那么當

時(8.14)式達到最大，稱為第一判別函數(shù)，而

是在約束條件

之下使得(8.14)式達到最大的解，稱為第二判別函數(shù)，如此下去

，是在約束條件

之下使得(8.14)式達到最大的解，稱為第s個判別函數(shù)．2024/7/31

18例8.2原油樣品數(shù)據(jù)的Fisher判別

2024/7/31

19表8-2三個沙巖層的原油樣品數(shù)據(jù)2024/7/31

20序號Gx1x2x3x4x5113.951.00.207.0612.19212.749.00.077.1412.23……………………………………613.943.00.076.2510.42712.735.00.005.119.00825.047.00.077.066.10923.432.00.205.824.69……………………………………1724.446.00.077.545.761823.030.00.005.1210.771936.313.00.504.248.272031.75.61.005.694.64……………………………………5535.034.00.704.216.505636.227.00.303.972.97例8.2(續(xù))原有樣品數(shù)據(jù)的Fisher判別>Z=predict(ld)>newG=Z$class>cbind(Species,newG,Z$post)#Z$post給出回判后驗概率>tab=table(newG,G)>tab

GnewG1231710

2091

30137對56個原始數(shù)據(jù)的回代判別中,有3個錯誤,誤判率為5.357%.2024/7/31

21例8.2(續(xù))原有樣品數(shù)據(jù)的Fisher判別>x=Z$x[,1];y=Z$x[,2];plot(Z$x,type="n")>points(x[1:7],y[1:7],pch=17)

#將總體G1用實心三角點標出>points(x[8:18],y[8:18],pch=25)

#將總體G2用空心倒三角點標出>points(x[19:56],y[19:56],pch=19)

#將總體G3用實心圓點標出>text(c(-4.2,-2,1),c(1.2,-1.8,0.1),labels=c("G1","G2","G3"))#在適當位置標出總體名稱>newdata=data.frame(x1=c(4.5,3.2,7.3),x2=c(33,43,22),x3=c(0.2,0.1,0.4),x4=c(6.5,6.8,4.6),x5=c(5.5,11.5,5.8))>(predict(ld,newdata))

#對3個新的原油樣品的類別進行判別$class#列出新樣品的判別分類[1]213Levels:1232024/7/31

222024/7/31

23圖8-456個原油樣品數(shù)據(jù)的Fisher回判結(jié)果分類圖8.3

Bayes判別Bayes判別法假定對研究對象已經(jīng)有一定的認識，這種認識可以用先驗概率來描述，當取得樣本后，就可以利用樣本來修正已有的先驗概率分布，得到后驗分布，再通過后驗分布進行各種統(tǒng)計推斷(事實上，F(xiàn)isher判別中也用到了后驗概率)．Bayes判別法屬于概率判別法．判別準則：(1)個體歸屬某類的概率最大；

(2)錯判總平均損失最小.2024/7/31

248.3.1

兩總體的Bayes判別設(shè)有兩個總體G1,G2.概率密度函數(shù)及先驗概率分別為

及判別準則是極小化平均誤判損失ECM：判別規(guī)則:樣品x被判入G1,若x∈R1,樣品x被判入G2,若x∈R2=Ω-R1對應(yīng)的樣本空間

的劃分為2024/7/31

258.3.2

多總體的Bayes判別設(shè)有k個總體G1,G2,…,

Gk

概率密度函數(shù)及先驗概率分別為

及2024/7/31

26并假設(shè)所有的錯判損失相同，對待判樣品x相應(yīng)的判別準則為（詳細討論參見教材）例8.3部分國家人文發(fā)展水平的Bayes

判別

聯(lián)合國2020年《人類發(fā)展報告》給出了人均國民收入x1(美元)、預(yù)期壽命x2(歲)和預(yù)期受教育年限x3(年)和平均受教育年限x4．現(xiàn)從2020年各國人文發(fā)展指數(shù)（簡稱HDI）排序中選取了極高發(fā)展水平、高發(fā)展水平和中等發(fā)展水平國家各六個作為三個已知分類樣品總體，另選四個國家日本、印度、中國和南非作為待判樣品(數(shù)據(jù)另見newdata8.3)，數(shù)據(jù)合并如表8－3所示．對此數(shù)據(jù)進行Bayes判別分析．2024/7/31

27表8-32020年部分國家人文發(fā)展水平和主要指標2024/7/31

28序號國家Gx1x2x3x41美國16382678.916.313.42德國15531481.317.014.23希臘13015582.217.910.64新加坡18815583.616.411.65意大利14277683.516.110.46韓國14304483.016.512.27古巴2862178.814.311.88伊朗21244776.714.810.39巴西21426375.915.48.010泰國21778177.215.07.911烏克蘭21321672.115.111.412印尼21145971.713.68.213尼泊爾3345770.812.85.014伊拉克31080170.611.37.315喀麥隆3358159.312.16.316巴基斯坦3500567.38.35.217緬甸3496167.110.75.018敘利亞3361372.78.95.119日本待判4293284.615.212.920印度待判668169.712.26.521中國待判1605776.914.08.122南非待判1212964.113.810.22024/7/31

29#例8.3三類不同人文發(fā)展水平國家樣品數(shù)據(jù)的Bayes判別分析>setwd("C:/data")#設(shè)定工作路徑>d8.3<-read.csv("exam8.3.csv",header=T)#將"exam8.3數(shù)據(jù)讀入>library(MASS)#加載MASS程序包，以便使用其中的lda()函數(shù)>ld=lda(G~x1+x2+x3+x4,data=d8.3,prior=c(1,1,1)/3);ld#作判別分析……>Z=predict(ld)#作回判預(yù)測>newG=Z$class#回判分類記作newG>cbind(G,newG,Z$post)#Bayes判別法把樣品判入后驗概率大的那一類GnewG1231119.999961e-013.946080e-062.483229e-20……6119.998596e-011.404463e-042.705233e-197228.702931e-059.999130e-012.167775e-09……18333.072589e-225.179978e-089.999999e-01>tab=table(newG,G);tab#列表比較>sum(diag(prop.table(tab)))#計算回判正確率#由程序輸出結(jié)果可見，三類水平共18個國家回判結(jié)果全部正確．2024/7/31

30>x=Z$x[,1];y=Z$x[,2]#取Z$x的兩列LD1和LD2構(gòu)成坐標面>plot(Z$x,type="n")>points(x[1:6],y[1:6],pch=24)#用空心上三角點標出總體G1樣品點>text(x[1:6],y[1:6],labels=exam8.3[1:6,2],adj=c(1.3,0.2),cex=0.7)

#標出第一類國家名>points(x[7:12],y[7:12],pch=1)#用空心圓點標出總體G2樣品點>text(x[7:12],y[7:12],labels=exam8.3[7:12,2],adj=c(-0.4,0.4),cex=0.7)

#標出第二類國家名>points(x[13:18],y[13:18],pch=25)#用空心倒三角點標出總體G3樣品點>text(x[13:18],y[13:18],labels=exam8.3[13:18,2],adj=c(-0.2,0.4),cex=0.7)

#標出第三類國家名2024/7/31

31>newdata8.3<-read.csv("newdata8.3.csv",header=T)

#將待判四國數(shù)據(jù)讀入>newdata=newdata8.3[1:4,4:7]>Z1=predict(ld,newdata);Z1#對4個待判國家所屬總體進行判別$class[1]1322Levels:123$posterior12319.996012e-010.00039879253.539070e-1821.976655e-140.00311291719.968871e-0133.996557e-060.99997069572.530774e-0542.943094e-100.82304747481.769525e-01>points(Z1$x[,1],Z1$x[,2],pch=19)#用實心圓點標出待判四國樣品點>newnames=c("日本","印度","中國","南非")>text(Z1$x[,1],Z1$x[,2],labels=newnames,adj=c(0.5,1.3),cex=0.7)

#標出四個國家名2024/7/31

32圖8-5三類人文發(fā)展水平國家及四個待判國家回判結(jié)果示意圖5.4二次判別2024/7/31

332024/7/31

34例8.4

iris3是鳶尾花數(shù)據(jù)iris的另一種形式，它將三種鳶尾花按品種分成三類羅列，每類50個數(shù)據(jù).2024/7/3135例8.4

(續(xù))現(xiàn)要從每類鳶尾花數(shù)據(jù)中各自無放回地隨機抽取40個數(shù)據(jù)，共120個數(shù)據(jù)，組成訓(xùn)練樣本集來建立二次判別函數(shù)，并利用它對剩下的30個樣本數(shù)據(jù)的類別進行二次判別．

2024/7/31

362024/7/31

378.5案例分析與

R實現(xiàn)

案例8.1(數(shù)據(jù)文件為case8.1)表8-4中列出了2020年頭7個月我國35個主要城市食品煙酒類城市居民消費價格指數(shù)CPI（上年同月=100）（%）．下面先利用前30個城市數(shù)據(jù)進行k=3的k均值聚類，再以這三個類為基礎(chǔ)，分別用Fisher判別法、Bayes判別法和距離判別法、對表中前30個城市進行回判，再對余下5個城市天津、海口、成都、昆明和烏魯木齊的屬類進行判別分析．2024/7/31

38表8-42020年頭7個月我國35個城市食品煙酒類城市居民消費價格指數(shù)2024/7/3139地區(qū)x1x2x3x4x5x6x7北京110.3109.8109.0106.9104.2105.3106.0石家莊113.0111.5110.6108.0107.0107.4109.5太原113.2113.8112.9108.6106.6107.1106.9呼和浩特107.9108.1106.8104.4102.0103.2103.7沈陽117.6116.0112.9109.2107.6106.1108.0大連113.0112.2111.0108.3106.7106.0107.6長春115.0113.9111.8110.4107.3104.8106.4…………………………………………西安110.8113.8109.1107.1105.2105.8105.6蘭州108.4110.3107.2106.4104.7104.5105.1西寧109.7111.7110.2107.8106.5107.2107.8銀川109.7110.2108.6107.3104.2103.4104.0天津111.8111.3110.3108.4106.4106.2107.3?？?14.5113.7113.5113.1107.1106.1106.2成都117.0121.5115.9113.8111.5113.2114.1昆明115.1118.2116.1112.9109.6107.8109.6烏魯木齊108.8106.0103.9103.4100.6102.9104.4（1）k均值聚類>setwd("C:/mdata")

#設(shè)定工作路徑>d8.1<-read.csv("case8.1.csv",header=T)

#將case8.1數(shù)據(jù)讀入>ca8.1=d8.1[,-1]

#d8.1的第一列為樣本名稱先去掉>rownames(ca8.1)=d8.1[,1]

#用d8.1的第一列為ca8.1的行重新命名>KM3<-kmeans(ca8.1[1:30,],3,nstart=15,algorithm="Hartigan-Wong")

#聚類數(shù)為3>sort(KM3$cluster)

#對分類結(jié)果進行排序并查看北京呼和浩特上海蘭州銀川石家莊太原大連長春杭州

111

1

1

2

2

2

22寧波福州南昌長沙西安西寧沈陽哈爾濱南京合肥

2222223333廈門濟南青島鄭州武漢廣州深圳南寧重慶貴陽

333333333

3>f8.1=sort(KM3$cluster);names(f8.1)>c8.1=ca8.1[names(f8.1),]

#對前30個城市的分類結(jié)果重新按類排序>G=rep(c(1,2,3),KM3[[7]])

#重復(fù)上述步驟時排序可能變化，用KM3[[7]]動態(tài)調(diào)整>cn8.1=cbind(G,c8.1);cn8.1#對新類指定類別號并展示

Gx1x2

x3x4x5x6x7北京1

110.3109.8109.0106.9104.2105.3106.0………………貴陽3115.5118.3115.5113.7110.2109.1111.92024/7/31

40（2）Fisher判別法>attach(cn8.1)#把數(shù)據(jù)變量名字放入內(nèi)存>library(MASS);>ld=lda(G~x1+x2+x3+x4+x5+x6+x7,data=cn8.1);ldCall:lda(G~x1+x2+x3+x4+x5+x6+x7,data=cn8.1)Priorprobabilitiesofgroups:

1

2

30.16666670.36666670.4666667Groupmeans:

x1x2x3x4x5

x6x71109.4400109.5600107.8400106.4400104.0800104.2800105.12002112.2091112.9182110.8636108.6000106.2909106.8182108.20913115.8929116.5643114.1786111.8214108.7714108.5929110.0286Coefficientsoflineardiscriminants:

LD1LD2x10.363916050.27863144……x70.11607995-0.11950152Proportionoftrace:

LD1

LD20.98280.01722024/7/31

41對原始數(shù)據(jù)進行回判分類，并與真實的分類進行對比>Z=predict(ld);newG=Z$class#進行回判，并記回判分類為newG>cbind(G,newG,Z$post)

#合并原分類、回判分類、回判后驗概率

GnewG1

2

3北京

1

19.329013e-016.709868e-02

1.517116e-09……銀川

119.998689e-011.311142e-048.634216e-14石家莊221.954933e-059.995079e-014.725858e-04……西寧

224.709723e-03

9.952887e-011.612069e-06沈陽3

35.199081e-148.354152e-04

9.991646e-01……貴陽

335.335377e-17

2.352881e-05

9.999765e-01>tab=table(G,newG);tab#原分類和新分類列表比較newGG12315002

0110300

14>sum(diag(prop.table(tab)))#計算判別符合率[1]12024/7/31

42三類城市的回判全部正確．再對5個待判城市（newdata）的屬類進行判別>newdata=ca8.1[31:35,]#選取待判別城市>predict(ld,newdata=newdata)#對5個待判城市的屬類進行判定$class[1]22331Levels:123$posterior

12

3天津

2.230797e-029.776865e-015.529849e-06海口

5.696717e-039.740976e-01

2.020564e-02成都

2.081615e-291.221230e-101.000000e+00昆明

3.670545e-168.264183e-059.999174e-01烏魯木齊

1.000000e+002.894338e-099.371758e-232024/7/31

43可以看出5個待判城市中:天津、海口被判入第2類；成都、昆明被判入第3類；烏魯木齊被判入第1類．（3）Bayes判別法>attach(cn8.1);library(MASS)#把數(shù)據(jù)變量名字放入內(nèi)存，并加載MASS包>ld=lda(G~x1+x2+x3+x4+x5+x6+x7,prior=c(1,1,1)/3,data=cn8.1);ldCall:lda(G~x1+x2+x3+x4+x5+x6+x7,data=cn8.1,prior=c(1,1,1)/3)Priorprobabilitiesofgroups:

1230.33333330.33333330.3333333Groupmeans:

x1x2x3x4x5

x6

x71109.4400109.5600107.8400106.4400104.0800104.2800105.12002112.2091112.9182110.8636108.6000106.2909106.8182108.20913115.8929116.5643114.1786111.8214108.7714108.5929110.0286Coefficientsoflineardiscriminants:

LD1

LD2x10.351594390.29402680……x70.12111268-0.11439794Proportionoftrace:

LD1

LD20.98560.01442024/7/31

44對原始數(shù)據(jù)進行回判分類，并與真實的分類進行對比>Z=predict(ld);newG=Z$class#進行回判，并記回判分類為newG>cbind(G,newG,Z$post)

#合并原分類、回判分類、回判后驗概率

GnewG1

2

3北京119.683419e-013.165806e-02

5.624109e-10……銀川

119.999404e-015.960161e-05

3.083869e-14石家莊

224.301187e-059.995856e-013.713463e-04……西寧221.030316e-029.896956e-011.259508e-06沈陽331.455411e-131.063014e-039.989370e-01……貴陽

331.493896e-162.994557e-059.999701e-01>tab=table(G,newG);tab#原分類和新分類列表比較newGG12315002

0110300

14>sum(diag(prop.table(tab)))#計算判別符合率[1]12024/7/31

45由輸出結(jié)果可見，三類城市的回判仍然全部正確，只是后驗概率的數(shù)值上有一定變化．再對5個待判城市（newdata）的屬類進行判別:>prenew=predict(ld,newdata=newdata);prenew#對五個待判樣本進行判定$class[1]22331Levels:1