《多元统计分析-基于R（第3版）》课件 第7、8章-聚类分析、判别分析

第7章聚类分析物以类聚，人以群分。引入统计和数学工具把研究对象进行合理分类十分重要．聚类分析（clusteranalysis）就是研究如何将由多个个体组成的研究对象按照个体之间的相似性进行合理分类的一种多元统计方法．聚类分析在经济、管理、医学、地质、生物等诸多领域内有广泛应用．本章介绍个体间的相似性度量及系统聚类、k均值聚类和EM聚类等常用的聚类方法．2024/7/312第7章聚类分析7.1相似性的度量各种距离：用于对样品进行聚类（Q型聚类）相关系数：用于对变量进行聚类（R型聚类）7.2

系统聚类法（一种基本的分层聚类法）7.3

k均值聚类法（一种常用快速聚类法）7.4

EM聚类法（一种迭代寻优聚类法）2024/7/3137.1相似性度量设是对p维总体的n次观测得到的n×p数据矩阵，n个行看成

n个样品的观测值，p个列看成p个变量的观测值．

xij是X的第i行第j列元素.42024/7/31Q型聚类对样品进行聚类,即根据样品间相互“靠近”的程度来进行聚类.通常用各种统计距离来描述这种靠近程度，距离小的样品聚为同一类，距离大的样品聚为不同的类．设与是两个样品点（可视为矩阵X的任意两行）.5

2024/7/31常用的几种距离欧式距离：绝对距离：切氏距离：6

2024/7/31明氏距离：马氏距离：兰氏距离：其中2024/7/31

7数据阵X=(xij)的列向量间相关系数

其中p为变量数，n为样品数.8

2024/7/31R型聚类对变量进行聚类，即根据变量间的“相似”度来进行聚类，可用变量（矩阵X的任意两列）间的相关系数或夹角余弦来描述这种相似性．数据阵X=(xij)的列向量间的夹角余弦：9

2024/7/317.2系统聚类法系统聚类法的步骤:(注意:对n个样品或者p个变量都可以进行系统聚类,只不过使用的相似性度量不同而已)(1)先将要分析的每个个体看成一类，共r类(对Q型聚类，r=n；对R型聚类，r=p);(2)找出最相似的两类，合并成一个新类，得r-1类;(3)在r-1类中，再找出最相似的两类合并，得r-2类;(4)以此类推，将所有的r类合并成一大类.10

2024/7/31设Gs和Gt为系统聚类过程中生成的任意两个子类，dij为Gs中第i样品与Gt中第j样品之间的距离，Dst为Gs与Gt之间的距离.这时,根据具体情况,可使用不同的距离度量来描述子类与子类之间的距离.(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个变量之间的相似系数，则Gs与Gt间的相似度通常用Rst来度量：

注意：也可以将变量间的相似系数cij转化成变量间的距离dij（例如

）来聚类.14

2024/7/313.R型聚类时子类之间相似程度的度量例7.1(数据文件为exam7.1)从湖南邓阜仙岩体采集了七块花岗岩样品，分别测得其五种化学成分（见表7-1），试用系统聚类的最小距离法和最大距离法对这七块花岗岩样品进行聚类．序号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随机采集的七块花岗岩样品的部分化学成分数据

解

首先采用最小距离法进行聚类，将七块花岗岩样品看成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系统聚类（最小距离法）R程序#example7.1系统聚类（假定数据存储目录为：c:/data）>setwd("C:/data")#设定工作路径>d7.1<-read.csv(“exam7.1.csv”,header=T)#将exam7.1数据读入到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为明氏距离参数k.>HC<-hclust(d,method="single")#采用最小距离法聚类#method为系统聚类方法,包括"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}.

函数cutree()可以将hclust()聚类结果按顺序标出样品所属类别：>cutree(HC1,k=3)

#指定分类个数k=3[1]11222337.3k均值聚类法1.基本思想

系统聚类法的每一步都要计算“类间距离”，计算量比较大.MacQueen(1967)提出了一种动态快速聚类方法—k均值聚类法(k-meanscluster).

其基本思想是：

根据给定的参数k,先把n个对象粗略的分为k类,然后按照某种最优原则(通常为一个准则函数)修改不合理的分类,直到准则函数收敛为止,就得到了一个最终的分类.

232024/7/31例7.2(数据文件为exam7.2)2024/7/3124

表7—5给出了2014-2020年金砖国家和七国集团总发电量数据（单位：太瓦时(10亿度)）．根据这些数据，采用

k均值聚类法进行聚类分析，k分别取4和5．2024/7/3125表7-5近七年来金砖国家和七国集团总发电量国家\年份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#金砖国家和七国集团近七年总发电量的k均值聚类分析．>setwd("C:/data")#设定工作路径>exam7.2<-read.csv("exam7.2.csv",header=T)#将exam7.2数据读入>d7.2=exam7.2[,-1]#exam7.2的第一列为国家名，不是数值先去掉>rownames(d7.2)=exam7.2[,1]#用exam7.2的第一列为d7.2的行重新命名>KM4<-kmeans(d7.2,4,nstart=20,algorithm=“Hartigan-Wong”)#聚类个数先取为4,初始随机集合个数取为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")#聚类个数取为5>sort(KM5$cluster)

#对分类结果进行排序并查看分类情况

中国

美国

印度俄罗斯日本

南非英国意大利

巴西加拿大

德国法国

1123344455552024/7/312727

2024/7/31聚为4类聚为5类中国中国，美国俄罗斯，印度，日本印度美国俄罗斯，日本南非，英国，意大利，巴西，加拿大，德国，法国，南非，英国，意大利巴西，加拿大，德国，法国分析比较：282024/7/31图7-3金砖国家和七国集团总发电量最小距离法聚类树状图7.4其他聚类函数2024/7/3129K-中心点聚类：pam()函数（需先加载clust软件包）密度聚类：dbscan()函数（需先下载加载fpc软件包）期望最大化EM聚类：(Expectation-Maximization)的Mclust()函数（需先下载并加载mclust软件包）等EM聚类函数Mclust()的使用格式：Mclust(data,G,modelNames,prior,control,...)其中data为待聚类数据集；G为预设类别数，默认值为1至9，由软件根据BIC值选择最优值；modelNames用于设定模型类别，也由函数自动选取最优值．2024/7/3130例7.3在R软件内置数据集中，有一个由地质学家于1978年8月至1979年8月在美国黄石公园旅游景点老忠实泉（OldFaithful）记录的间歇喷泉喷发数据集，名为faithful．数据集有272行，2列，两列数据各为泉水喷发持续时间（eruptions）和喷发间隔时间（waiting），时间单位均为分钟．

以下用两种方法进行EM聚类：（1）用函数Mclust直接进行EM聚类：>library(mclust)#加载mclust软件包>EM1<-Mclust(faithful)#直接做EM聚类

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

#绘制聚类的概率分布图

（全部272个原始数据被聚为三类,详细分析过程参见教材）31

2024/7/31图7-4

对faithful数据直接进行EM聚类结果图2024/7/3132（2）首先在faithful数据分布范围内随机生成728个均匀分布随机数，将它们与原来的faithful数据混合得到大小为1000的混合样本数据，再用函数Mclust对混合样本作EM聚类分析．>nNoise<-728#设定均匀分布噪声数据个数>set.seed(9)#设置随机数种子>Noise<-apply(faithful,2,function(x)runif(nNoise,min=min(x)-0.1,max=max(x)+0.1))#在faithful数据分布范围内生成nNoise=728行,2列的均匀分布噪声数据>data<-rbind(faithful,Noise)#按行合并faithful和Noise，得到1000个混合数据样本>plot(faithful)

#绘制喷发-间隔数据散点图>points(Noise,pch=16,cex=0.5)

#在散点图中汇入噪声数据点>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数据与均匀分布数据的混合样本进行EM聚类的结果图2024/7/3134

从程序输出结果和图7-5易见，全部1000个原始数据被聚为2类，样本大小分别为83(蓝圆点)和167(红空心方块点)，其余750个点均被视为噪声点．两个类的均值分别为（2.1209,54.1188）和（4.3312,79.9068），同样也输出了两类数据对应的协方差矩阵以及似然函数值和BIC值等．

从图7-5还可以看出，对混合样本的EM聚类基本没有受到均匀分布噪声的影响，将faithful数据聚为2类.这与（1）中直接聚类法聚成3类的结果有所区别，所以不同的聚类方法会产生不同的聚类结果．2024/7/3135

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

37判别分析是在已知样品所有可能分类的前提下，将给定的新样品按照某种分类准则判入其中某个类中的一种多元统计方法．例如：根据患者的各项检查指标来判断该病人属于哪类病症；根据某地气象的记录资料来判别(预报)未来几天的天气状况；根据某地相关经济指标判断该地区属于哪一种经济类型地区．2024/7/318.1距离判别距离判别简介

、

两个及多个总体的距离判别8.2

Fisher判别两个及多个总体的Fisher判别8.3

Bayes判别两个及多个总体的Bayes判别8.4二次判别8.5案例分析与R实现本章主要内容：388.1距离判别思想：根据一个样品与各个类别距离的远近对

其所属类别进行判定．这里的距离通常使用马氏距离：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.试据此数据做距离判别分析．今若测得两个患者的两项指标分别为(90,160)和(85,155)，对他们如何进行判断？冠心病人正常人组别x1x2组别x1x2174200294172110014421001181110150270152170274280172196212280190180158270142180172280107110014028012411002302801941100220278152190239270190111015528010411001552809419614028413211002302701402024/7/31

44>setwd("C:/data")#设定工作路径>d8.1<-read.csv("exam8.1.csv",header=T)#将exam8.1数据读入>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#计算两个总体样的协方差矩阵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年在植物分类研究中提出了该判别法．主要思想是通过将多维数据投影到一维直线上，使得同一类别(总体)中的数据在该直线上尽量靠拢，不同类别(总体)的数据尽可能分开，然后再利用前面的距离判别法来建立判别准则．这种投影直线（判别函数）可能有一条或多条．常用方法：线性判别法；非线性判别法；典型判别法等．这里主要介绍线性判别法．2024/7/31

485.2.1两总体的Fisher判别2024/7/31

49

两总体Fisher判别的思想是将二维空间中的点投影到一条适当的直线

y

上，使得两个总体G1和

G2中的点在直线y上的投影点尽可能分开，而同一总体在y上的投影点尽可能靠拢，再利用上面的距离判别法来建立判别准则．

示意图：总体G1（大圆点）和

G2

（小圆点）中的点在直线y上的投影点满足上面要求．再根据新样品点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判别多总体情形，通常要选取多条投影直线，即选取多个判别函数

来进行判别．设有k个总体G1,G2,···,

Gk，它们有共同的协方差阵Σ，均值分别为μ1,μ2

,···,μk.令问题为选择a,使的投影变差比值达到最大:2024/7/31

528.2.2

多总体Fisher判别定理8.2设为的s个非零特征根，为相应的特征向量且满足

，那么当

时(8.14)式达到最大，称为第一判别函数，而

是在约束条件

之下使得(8.14)式达到最大的解，称为第二判别函数，如此下去

，是在约束条件

之下使得(8.14)式达到最大的解，称为第s个判别函数．2024/7/31

53例8.2原油样品数据的Fisher判别

2024/7/31

54表8-2三个沙岩层的原油样品数据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(续)原有样品数据的Fisher判别>Z=predict(ld)>newG=Z$class>cbind(Species,newG,Z$post)#Z$post给出回判后验概率>tab=table(newG,G)>tab

GnewG1231710

2091

30137对56个原始数据的回代判别中,有3个错误,误判率为5.357%.2024/7/31

56例8.2(续)原有样品数据的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个原油样品数据的Fisher回判结果分类图8.3

Bayes判别Bayes判别法假定对研究对象已经有一定的认识，这种认识可以用先验概率来描述，当取得样本后，就可以利用样本来修正已有的先验概率分布，得到后验分布，再通过后验分布进行各种统计推断(事实上，Fisher判别中也用到了后验概率)．Bayes判别法属于概率判别法．判别准则：(1)个体归属某类的概率最大；

(2)错判总平均损失最小.2024/7/31

598.3.1

两总体的Bayes判别设有两个总体G1,G2.概率密度函数及先验概率分别为

及判别准则是极小化平均误判损失ECM：判别规则:样品x被判入G1,若x∈R1,样品x被判入G2,若x∈R2=Ω-R1对应的样本空间

的划分为2024/7/31

608.3.2

多总体的Bayes判别设有k个总体G1,G2,…,

Gk

概率密度函数及先验概率分别为

及2024/7/31

61并假设所有的错判损失相同，对待判样品x相应的判别准则为（详细讨论参见教材）例8.3部分国家人文发展水平的Bayes

判别

联合国2020年《人类发展报告》给出了人均国民收入x1(美元)、预期寿命x2(岁)和预期受教育年限x3(年)和平均受教育年限x4．现从2020年各国人文发展指数（简称HDI）排序中选取了极高发展水平、高发展水平和中等发展水平国家各六个作为三个已知分类样品总体，另选四个国家日本、印度、中国和南非作为待判样品(数据另见newdata8.3)，数据合并如表8－3所示．对此数据进行Bayes判别分析．2024/7/31

62表8-32020年部分国家人文发展水平和主要指标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三类不同人文发展水平国家样品数据的Bayes判别分析>setwd("C:/data")#设定工作路径>d8.3<-read.csv("exam8.3.csv",header=T)#将"exam8.3数据读入>library(MASS)#加载MASS程序包，以便使用其中的lda()函数>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)

#将待判四国数据读入>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三类人文发展水平国家及四个待判国家回判结果示意图5.4二次判别2024/7/31

682024/7/31

69例8.4

iris3是鸢尾花数据iris的另一种形式，它将三种鸢尾花按品种分成三类罗列，每类50个数据.2024/7/3170例8.4

(续)现要从每类鸢尾花数据中各自无放回地随机抽取40个数据，共120个数据，组成训练样本集来建立二次判别函数，并利用它对剩下的30个样本数据的类别进行二次判别．

2024/7/31

712024/7/31

728.5案例分析与

R实现

案例8.1(数据文件为case8.1)表8-4中列出了2020年头7个月我国35个主要城市食品烟酒类城市居民消费价格指数CPI（上年同月=100）（%）．下面先利用前30个城市数据进行k=3的k均值聚类，再以这三个类为基础，分别用Fisher判别法、Bayes判别法和距离判别法、对表中前30个城市进行回判，再对余下5个城市天津、海口、成都、昆明和乌鲁木齐的属类进行判别分析．2024/7/31

73表8-42020年头7个月我国35个城市食品烟酒类城市居民消费价格指数2024/7/3174地区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海口114.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数据读入>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")

#聚类数为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]]动态调整>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)#把数据变量名字放入内存>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对原始数据进行回判分类，并与真实的分类进行对比>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海口

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个待判城市中:天津、海口被判入第2类；成都、昆明被判入第3类；乌鲁木齐被判入第1类．（3）Bayes判别法>attach(cn8.1);library(MASS)#把数据变量名字放入内存，并加载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对原始数据进行回判分类，并与真实的分类进行对比>Z=predict(ld);newG=Z$class#进行回判，并记回判分类为newG>cbind(G,newG,Z$post)

#合并原分类、回判分类、回判后验概率

GnewG1

2

3北京119.683419e-013.165806e-02

5.624109e-10