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

主编：费宇中国人民大学出版社第8章判别分析2024/7/31

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

、

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

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

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

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

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

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年在植物分类研究中提出了该判别法．主要思想是通过将多维数据投影到一维直线上，使得同一类别(总体)中的数据在该直线上尽量靠拢，不同类别(总体)的数据尽可能分开，然后再利用前面的距离判别法来建立判别准则．这种投影直线（判别函数）可能有一条或多条．常用方法：线性判别法；非线性判别法；典型判别法等．这里主要介绍线性判别法．2024/7/31

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

14

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

y

上，使得两个总体G1和

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

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

G2

（小圆点）中的点在直线y上的投影点满足上面要求．再根据新样品点x在

y上的投影点

y=aTx靠近μ1y（在μy右侧）还是靠近μ2y（在μy左侧）来判断新样品点x的归属．2024/7/31

15图8-3两总体Fisher判别示意图两总体的Fisher判别准则：2024/7/31

16

8.2.2

多总体Fisher判别多总体情形，通常要选取多条投影直线，即选取多个判别函数

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

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

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

178.2.2

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

，那么当

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

是在约束条件

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

，是在约束条件

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

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

2024/7/31

19表8-2三个沙岩层的原油样品数据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(续)原有样品数据的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

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

222024/7/31

23图8-456个原油样品数据的Fisher回判结果分类图8.3

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

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

248.3.1

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

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

的划分为2024/7/31

258.3.2

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

Gk

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

及2024/7/31

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

判别

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

27表8-32020年部分国家人文发展水平和主要指标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三类不同人文发展水平国家样品数据的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

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

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

32图8-5三类人文发展水平国家及四个待判国家回判结果示意图5.4二次判别2024/7/31

332024/7/31

34例8.4

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

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

2024/7/31

362024/7/31

378.5案例分析与

R实现

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

38表8-42020年头7个月我国35个城市食品烟酒类城市居民消费价格指数2024/7/3139地区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

40（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

41对原始数据进行回判分类，并与真实的分类进行对比>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)#把数据变量名字放入内存，并加载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对原始数据进行回判分类，并与真实的分类进行对比>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由输出结果可见，三类城市的回判仍然全部正确，只是后验概率的数值上有一定变化．再对5个待判城市（newdata）的属类进行判别:>prenew=predict(ld,newdata=newdata);prenew#对五个待判样本进行判定$class[1]22331Levels:1