《多元統(tǒng)計分析-基于R（第3版）》課件 第7、8章-聚類分析、判別分析

第7章聚類分析物以類聚，人以群分。引入統(tǒng)計和數(shù)學工具把研究對象進行合理分類十分重要．聚類分析（clusteranalysis）就是研究如何將由多個個體組成的研究對象按照個體之間的相似性進行合理分類的一種多元統(tǒng)計方法．聚類分析在經(jīng)濟、管理、醫(yī)學、地質、生物等諸多領域內有廣泛應用．本章介紹個體間的相似性度量及系統(tǒng)聚類、k均值聚類和EM聚類等常用的聚類方法．2024/7/312第7章聚類分析7.1相似性的度量各種距離：用于對樣品進行聚類（Q型聚類）相關系數(shù)：用于對變量進行聚類（R型聚類）7.2

系統(tǒng)聚類法（一種基本的分層聚類法）7.3

k均值聚類法（一種常用快速聚類法）7.4

EM聚類法（一種迭代尋優(yōu)聚類法）2024/7/3137.1相似性度量設是對p維總體的n次觀測得到的n×p數(shù)據(jù)矩陣，n個行看成

n個樣品的觀測值，p個列看成p個變量的觀測值．

xij是X的第i行第j列元素.42024/7/31Q型聚類對樣品進行聚類,即根據(jù)樣品間相互“靠近”的程度來進行聚類.通常用各種統(tǒng)計距離來描述這種靠近程度，距離小的樣品聚為同一類，距離大的樣品聚為不同的類．設與是兩個樣品點（可視為矩陣X的任意兩行）.5

2024/7/31常用的幾種距離歐式距離：絕對距離：切氏距離：6

2024/7/31明氏距離：馬氏距離：蘭氏距離：其中2024/7/31

7數(shù)據(jù)陣X=(xij)的列向量間相關系數(shù)

其中p為變量數(shù)，n為樣品數(shù).8

2024/7/31R型聚類對變量進行聚類，即根據(jù)變量間的“相似”度來進行聚類，可用變量（矩陣X的任意兩列）間的相關系數(shù)或夾角余弦來描述這種相似性．數(shù)據(jù)陣X=(xij)的列向量間的夾角余弦：9

2024/7/317.2系統(tǒng)聚類法系統(tǒng)聚類法的步驟:(注意:對n個樣品或者p個變量都可以進行系統(tǒng)聚類,只不過使用的相似性度量不同而已)(1)先將要分析的每個個體看成一類，共r類(對Q型聚類，r=n；對R型聚類，r=p);(2)找出最相似的兩類，合并成一個新類，得r-1類;(3)在r-1類中，再找出最相似的兩類合并，得r-2類;(4)以此類推，將所有的r類合并成一大類.10

2024/7/31設Gs和Gt為系統(tǒng)聚類過程中生成的任意兩個子類，dij為Gs中第i樣品與Gt中第j樣品之間的距離，Dst為Gs與Gt之間的距離.這時,根據(jù)具體情況,可使用不同的距離度量來描述子類與子類之間的距離.(1)最小距離法：11

2024/7/312.Q型聚類過程中子類之間相似程度的度量(2)最大距離法：(3)中間距離法：

其中(4)重心距離法：其中,和分別表示Gs和Gt的重心.12

2024/7/31(5)類平均距離法：(6)離差平方和法(Ward法)：

其中,13

2024/7/31設Gs和Gt為兩個子類,用rij表示Gs中第i個變量與Gt中第j個變量之間的相似系數(shù)，則Gs與Gt間的相似度通常用Rst來度量：

注意：也可以將變量間的相似系數(shù)cij轉化成變量間的距離dij（例如

）來聚類.14

2024/7/313.R型聚類時子類之間相似程度的度量例7.1(數(shù)據(jù)文件為exam7.1)從湖南鄧阜仙巖體采集了七塊花崗巖樣品，分別測得其五種化學成分（見表7-1），試用系統(tǒng)聚類的最小距離法和最大距離法對這七塊花崗巖樣品進行聚類．序號SiO2TiO2FeOCaOK2O175.200.141.860.915.21275.150.162.110.744.93372.190.131.520.694.65472.350.131.370.834.87572.740.101.410.724.99673.290.0331.070.173.15773.720.0330.770.282.782024/7/3115表7-1隨機采集的七塊花崗巖樣品的部分化學成分數(shù)據(jù)

解

首先采用最小距離法進行聚類，將七塊花崗巖樣品看成7個基本類，它們之間的距離（R計算程序見后）如表7—2所示．樣品123456710.000

20.4160.000

33.0883.0320.000

42.9132.8980.3410.000

52.5182.5110.6570.4260.000

63.0122.8371.9862.0922.0270.000

73.1132.9482.5652.6302.5400.6510.0002024/7/3116表7-2七塊花崗巖樣品之間的歐氏距離表4-2七塊花崗巖樣品按最小距離法的合并順序172024/7/31合并

次序要合并

的類合并后的新類最小距離法合并距離(歐氏距離)10.34120.41630.42640.65151.98662.511表7－3按最小距離法的合并順序及合并距離表4-3

七塊花崗巖樣品按最大距離法的合并順序18

2024/7/31合并

次序要合并

的類合并后的新類最大距離法合并距離(歐氏距離)10.34120.41630.65140.65752.63063.113從表7-3和表7-4可見：兩種聚類方法合并次序第1、2、5、6步相同，但第3、4兩步剛好相反．表7－4按最大距離法的合并順序及合并距離19

2024/7/31系統(tǒng)聚類（最小距離法）R程序#example7.1系統(tǒng)聚類（假定數(shù)據(jù)存儲目錄為：c:/data）>setwd("C:/data")#設定工作路徑>d7.1<-read.csv(“exam7.1.csv”,header=T)#將exam7.1數(shù)據(jù)讀入到d7.1中>d<-dist(d7.1,method="euclidean",diag=T,upper=F,p=2)#采用歐氏距離計算距離矩陣d,method為距離計算方法,包括“euclidean”（歐氏距離）,“manhattan”（絕對距離）,“maximum”（切氏距離）,“minkowski”（明氏距離）,"canberra"（蘭氏距離）等;diag為是否輸出對角線上的值,upper為是否輸出d的上三角部分的值,p為明氏距離參數(shù)k.>HC<-hclust(d,method="single")#采用最小距離法聚類#method為系統(tǒng)聚類方法,包括"single"（最小距離法）,"complete"（最大距離法）,"average"（類平均法）,"median"（中間距離法）,"centroid"（重心法）,"ward"（Ward法）等>plot(HC,hang=-1)#繪制最小距離法聚類樹狀圖（圖7-1）.

#當hang取負值時，從底部對齊開始繪制聚類樹狀圖>abline(h=c(2.2,1),lty=3)#在圖7-1中分別畫合并距離為2.2和1的兩條水平虛線，用來幫助我們確定分類202024/7/31圖7-1

最小距離法樹狀圖圖7-2

最大距離法樹狀圖從表7-1和圖7-1可以看出：若取合并距離為2.2（上虛線）,則7塊樣品可以分為兩類:第一類為{1,2},第二類為{3,4,5,6,7}.取合并距離為1（下虛線）,則7塊樣品可以分為三類:第一類為{1,2},第二類為{3,4,5},

第三類為{6,7}.2024/7/3121>HC1<-hclust(d,method="complete")#采用最大距離法聚類>x11()

#另開一個繪圖窗口>plot(HC1,hang=-1)#繪制最大距離法聚類樹狀圖（圖7-2）>rect.hclust(HC1,k=3,border="red")#用紅色矩形框出3個分類222024/7/31最大距離法聚類的R程序由此可知：第一類為{1,2},第二類為{3,4,5},第三類為{6,7}.

函數(shù)cutree()可以將hclust()聚類結果按順序標出樣品所屬類別：>cutree(HC1,k=3)

#指定分類個數(shù)k=3[1]11222337.3k均值聚類法1.基本思想

系統(tǒng)聚類法的每一步都要計算“類間距離”，計算量比較大.MacQueen(1967)提出了一種動態(tài)快速聚類方法—k均值聚類法(k-meanscluster).

其基本思想是：

根據(jù)給定的參數(shù)k,先把n個對象粗略的分為k類,然后按照某種最優(yōu)原則(通常為一個準則函數(shù))修改不合理的分類,直到準則函數(shù)收斂為止,就得到了一個最終的分類.

232024/7/31例7.2(數(shù)據(jù)文件為exam7.2)2024/7/3124

表7—5給出了2014-2020年金磚國家和七國集團總發(fā)電量數(shù)據(jù)（單位：太瓦時(10億度)）．根據(jù)這些數(shù)據(jù)，采用

k均值聚類法進行聚類分析，k分別取4和5．2024/7/3125表7-5近七年來金磚國家和七國集團總發(fā)電量國家\年份2014201520162017201820192020中國5794.55814.66133.26604.47166.17503.47779.1俄羅斯1064.21067.51091.01091.21109.21118.11085.4印度1262.21317.31401.71471.31579.21603.71560.9巴西590.5581.2578.9589.3601.4626.3620.1南非254.8250.1252.7255.1256.3252.6239.5美國4363.34348.74347.94302.54461.64411.24286.6加拿大647.6659.3663.7660.1655.8648.7643.9德國627.8647.6649.7652.9642.9609.4571.9法國564.9571.8556.2554.0574.0562.8524.9英國338.1338.9339.2338.2332.8324.8312.8意大利279.8283.0289.8295.8289.7293.9282.7日本1062.71030.11035.11042.11053.21030.31004.8#金磚國家和七國集團近七年總發(fā)電量的k均值聚類分析．>setwd("C:/data")#設定工作路徑>exam7.2<-read.csv("exam7.2.csv",header=T)#將exam7.2數(shù)據(jù)讀入>d7.2=exam7.2[,-1]#exam7.2的第一列為國家名，不是數(shù)值先去掉>rownames(d7.2)=exam7.2[,1]#用exam7.2的第一列為d7.2的行重新命名>KM4<-kmeans(d7.2,4,nstart=20,algorithm=“Hartigan-Wong”)#聚類個數(shù)先取為4,初始隨機集合個數(shù)取為20,算法為"Hartigan-Wong">sort(KM4$cluster)#對分類結果進行排序并查看中國俄羅斯印度日本美國巴西南非加拿大德國法國英國意大利

1222

344444442024/7/312626

2024/7/31

k均值聚類R程序>KM5<-kmeans(d7.2,5,nstart=10,algorithm="Hartigan-Wong")#聚類個數(shù)取為5>sort(KM5$cluster)

#對分類結果進行排序并查看分類情況

中國

美國

印度俄羅斯日本

南非英國意大利

巴西加拿大

德國法國

1123344455552024/7/312727

2024/7/31聚為4類聚為5類中國中國，美國俄羅斯，印度，日本印度美國俄羅斯，日本南非，英國，意大利，巴西，加拿大，德國，法國，南非，英國，意大利巴西，加拿大，德國，法國分析比較：282024/7/31圖7-3金磚國家和七國集團總發(fā)電量最小距離法聚類樹狀圖7.4其他聚類函數(shù)2024/7/3129K-中心點聚類：pam()函數(shù)（需先加載clust軟件包）密度聚類：dbscan()函數(shù)（需先下載加載fpc軟件包）期望最大化EM聚類：(Expectation-Maximization)的Mclust()函數(shù)（需先下載并加載mclust軟件包）等EM聚類函數(shù)Mclust()的使用格式：Mclust(data,G,modelNames,prior,control,...)其中data為待聚類數(shù)據(jù)集；G為預設類別數(shù)，默認值為1至9，由軟件根據(jù)BIC值選擇最優(yōu)值；modelNames用于設定模型類別，也由函數(shù)自動選取最優(yōu)值．2024/7/3130例7.3在R軟件內置數(shù)據(jù)集中，有一個由地質學家于1978年8月至1979年8月在美國黃石公園旅游景點老忠實泉（OldFaithful）記錄的間歇噴泉噴發(fā)數(shù)據(jù)集，名為faithful．數(shù)據(jù)集有272行，2列，兩列數(shù)據(jù)各為泉水噴發(fā)持續(xù)時間（eruptions）和噴發(fā)間隔時間（waiting），時間單位均為分鐘．

以下用兩種方法進行EM聚類：（1）用函數(shù)Mclust直接進行EM聚類：>library(mclust)#加載mclust軟件包>EM1<-Mclust(faithful)#直接做EM聚類

>summary(EM1,parameter=TRUE)#查看聚類結果>plot(EM1,what=“classification”)

#繪制聚類的概率分布圖

（全部272個原始數(shù)據(jù)被聚為三類,詳細分析過程參見教材）31

2024/7/31圖7-4

對faithful數(shù)據(jù)直接進行EM聚類結果圖2024/7/3132（2）首先在faithful數(shù)據(jù)分布范圍內隨機生成728個均勻分布隨機數(shù)，將它們與原來的faithful數(shù)據(jù)混合得到大小為1000的混合樣本數(shù)據(jù)，再用函數(shù)Mclust對混合樣本作EM聚類分析．>nNoise<-728#設定均勻分布噪聲數(shù)據(jù)個數(shù)>set.seed(9)#設置隨機數(shù)種子>Noise<-apply(faithful,2,function(x)runif(nNoise,min=min(x)-0.1,max=max(x)+0.1))#在faithful數(shù)據(jù)分布范圍內生成nNoise=728行,2列的均勻分布噪聲數(shù)據(jù)>data<-rbind(faithful,Noise)#按行合并faithful和Noise，得到1000個混合數(shù)據(jù)樣本>plot(faithful)

#繪制噴發(fā)-間隔數(shù)據(jù)散點圖>points(Noise,pch=16,cex=0.5)

#在散點圖中匯入噪聲數(shù)據(jù)點>NoiseInit<-sample(c(TRUE,FALSE),size=nrow(faithful)+nNoise,replace=TRUE,prob=c(3,1)/4)>EM2<-Mclust(data,initialization=list(noise=NoiseInit))

#EM聚類>summary(EM2,parameter=TRUE)

#查看模型建模結果>plot(EM2,what="classification")

#繪制聚類結果的概率分布圖33

2024/7/31圖7-5

對faithful數(shù)據(jù)與均勻分布數(shù)據(jù)的混合樣本進行EM聚類的結果圖2024/7/3134

從程序輸出結果和圖7-5易見，全部1000個原始數(shù)據(jù)被聚為2類，樣本大小分別為83(藍圓點)和167(紅空心方塊點)，其余750個點均被視為噪聲點．兩個類的均值分別為（2.1209,54.1188）和（4.3312,79.9068），同樣也輸出了兩類數(shù)據(jù)對應的協(xié)方差矩陣以及似然函數(shù)值和BIC值等．

從圖7-5還可以看出，對混合樣本的EM聚類基本沒有受到均勻分布噪聲的影響，將faithful數(shù)據(jù)聚為2類.這與（1）中直接聚類法聚成3類的結果有所區(qū)別，所以不同的聚類方法會產生不同的聚類結果．2024/7/3135

第8章判別分析2024/7/31

37判別分析是在已知樣品所有可能分類的前提下，將給定的新樣品按照某種分類準則判入其中某個類中的一種多元統(tǒng)計方法．例如：根據(jù)患者的各項檢查指標來判斷該病人屬于哪類病癥；根據(jù)某地氣象的記錄資料來判別(預報)未來幾天的天氣狀況；根據(jù)某地相關經(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)本章主要內容：388.1距離判別思想：根據(jù)一個樣品與各個類別距離的遠近對

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

2024/7/3140

2024/7/3141

2024/7/3142

2024/7/31例8.1

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

43

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

44>setwd("C:/data")#設定工作路徑>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

45>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

46圖8－2

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

2024/7/31

478.2

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

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

49

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

y

上，使得兩個總體G1和

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

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

G2

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

y上的投影點

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

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

51

8.2.2

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

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

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

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

528.2.2

多總體Fisher判別定理8.2設為的s個非零特征根，為相應的特征向量且滿足

，那么當

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

是在約束條件

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

，是在約束條件

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

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

2024/7/31

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

55序號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

56例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

572024/7/31

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

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

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

598.3.1

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

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

的劃分為2024/7/31

608.3.2

多總體的Bayes判別設有k個總體G1,G2,…,

Gk

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

及2024/7/31

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

判別

聯(lián)合國2020年《人類發(fā)展報告》給出了人均國民收入x1(美元)、預期壽命x2(歲)和預期受教育年限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

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

63序號國家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

64#例8.3三類不同人文發(fā)展水平國家樣品數(shù)據(jù)的Bayes判別分析>setwd("C:/data")#設定工作路徑>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)#作回判預測>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)))#計算回判正確率#由程序輸出結果可見，三類水平共18個國家回判結果全部正確．2024/7/31

65>x=Z$x[,1];y=Z$x[,2]#取Z$x的兩列LD1和LD2構成坐標面>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

66>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

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

682024/7/31

69例8.4

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

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

2024/7/31

712024/7/31

728.5案例分析與

R實現(xiàn)

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

73表8-42020年頭7個月我國35個城市食品煙酒類城市居民消費價格指數(shù)2024/7/3174地區(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")

#設定工作路徑>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)

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

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個城市的分類結果重新按類排序>G=rep(c(1,2,3),KM3[[7]])

#重復上述步驟時排序可能變化，用KM3[[7]]動態(tài)調整>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

75（2）Fisher判別法>attach(cn8.1)#把數(shù)據(jù)變量名字放入內存>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

76對原始數(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

77三類城市的回判全部正確．再對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?？?/p>

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

78可以看出5個待判城市中:天津、?？诒慌腥氲?類；成都、昆明被判入第3類；烏魯木齊被判入第1類．（3）Bayes判別法>attach(cn8.1);library(MASS)#把數(shù)據(jù)變量名字放入內存，并加載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

79對原始數(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