多元统计判别分析试验报告

1、判别分析（设计性实验）（Discriminantanalysis）实验原理：判别分析是判别样品所属类型的一种统计方法。判别分析是在已知研究对象分成若干类型（或组别）并已取得各种类型的一批已知样品的观测数目，在此基础上根据某些准则建立判别式，然后对未知类型的样品进行判别分类。本实验要求学生应用距离判别准则（即，对任给的一次观测，若它与第i类的重心距离最近，就认为它来自第i类），对两总体和多总体情形下分别进行判别分析。实验中需注意协方差矩阵相等时，选取线性判别函数；协方差矩阵不相等时，应选取二次判别函数。实验题目一：为了检测潜在的血友病A携带者，下表中给出了两组数据：（t11a8）非携带者（nj被

2、迫携带者（IL）GroupXiX2GroupXiX2i-0.0056-0.i6572-0.34780.ii5ii-0.i698-0.i5852-0.36i8-0.2008i-0.3469-0.i8792-0.4986-0.086i-0.08940.00642-0.50i5-0.2984i-0.i6790.07i32-0.i3260.0097i-0.08360.0i062-0.69ii-0.339i-0.i979-0.00052-0.36080.i237i-0.07620.03922-0.4535-0.i682i-0.i9i3-0.2i232-0.3479-0.i72ii-0.i092-0.ii

3、92-0.35390.0722i-0.5268-0.47732-0.47i9-0.i079i-0.08420.02482-0.36i-0.0399i-0.0225-0.0582-0.32260.i67i0.00840.07822-0.43i9-0.0687i-0.i827-0.ii382-0.2734-0.002i0.i2370.2i42-0.55730.05481-0.4702-0.30992-0.3755-0.18651-0.1519-0.06862-0.495-0.015310.0006-0.11532-0.5107-0.24831-0.2015-0.04982-0.16520.2132

4、1-0.1932-0.22932-0.2447-0.040710.15070.09332-0.4232-0.09981-0.1259-0.06692-0.23750.28761-0.1551-0.12322-0.22050.00461-0.1952-0.10072-0.2154-0.021910.02910.04422-0.34470.00971-0.228-0.1712-0.254-0.05731-0.0997-0.07332-0.3778-0.26821-0.1972-0.06072-0.4046-0.11621-0.0867-0.0562-0.06390.15692-0.3351-0.1

5、3682-0.01490.15392-0.03120.142-0.174-0.07762-0.14160.16422-0.15080.11372-0.09640.05312-0.26420.08672-0.02340.08042-0.33520.08752-0.18780.2512-0.17440.18922-0.4055-0.24182-0.24440.16142-0.47840.0282其中x1=log10(AHFactivity),x2=log10(AHFantigen)。卜表给出了五个新的观测，试对这些观测判别归类；(t11b8)观测%X21-.112-0.2792-.059-0.06

6、83.0640.0124-.043-0.0525-.050-0.098实验要求：(1)分别检验两组数据是否大致满足二元正态性；(2)分别计算两组数据的协方差矩阵，是否可以认为两者近似相等？(3)对训练样本和新观测合并作散点图，不同的类用不同颜色标识；(4)用Ida函数做判别分析，即在协方差矩阵相等的情形下作判别分析；(5)用qda函数做判别分析，即在协方差矩阵不相等的情形下作判别分析;(6)比较方法(4)和方法(5)的误判率。实验题目二:某商学研究生院的招生官员利用指标一一大学期间平均成绩GPA和研究生管理能力#试GMAT出了三类申请者的的成绩，将申请者分为三类：接受，不接受GPA与GMAT成

7、绩：(t11a6)4土之)彳寸Ato卜表中给GPA(xl)GMAT(x2)接受GPAGMAT(x1)(x2)不接受GPA(x1)GMAT(x2)待定2.9659612.5444622.8649433.1447312.4342522.8549633.2248212.247423.1441933.2952712.3653123.2837133.6950512.5754222.8944733.4669312.3540623.1531333.0362612.5141223.540233.1966312.5145822.8948533.6344712.3639922.844433.5958812.364

8、8223.1341633.356312.6642023.0147133.455312.6841422.7949033.557212.4853322.8943133.7859112.4650922.9144633.4469212.6350422.7554633.4852812.4433622.7346733.4755212.1340823.1246333.3552012.4146923.0844033.3954312.5553823.0341933.2852312.315052350933.2153012.4148923.0343833.5856412.1941123.0539933.33565

9、12.3532122.8548333.443112.639423.0145333.3860512.5552823.0341433.2666412.7239923.0444633.660912.8538123.3755912.938423.852113.7664613.244671实验要求：(1)对上表中的数据作散点图，不同的类用不同的颜色标识；(2)用Ida函数做判别分析，即在协方差矩阵相等的情形下作判别分析；(3)用qda函数做判别分析，即在协方差矩阵不相等的情形下作判别分析；(4)比较方法(2)和方法(3)的误判率；(5)现有一新申请者的GPA为3.21,GMAT成绩为497。请将该观测在

10、(1)的散点图中标出，并分别用方法(2)和方法(3)将其归类？你认为哪一种方法更合适？(6)观察(1)的散点图中第三类的观测点有无异常值？若有，将该异常值剔除后再对新申请者判别归类，结果有无变化？实验题目一分析报告：输出结果及分析：(1)分别检验两组数据是否大致满足二元正态性；>data5<-read.csv("data5.csv",header=T)#导入数据>data5b<-read.csv("data5b.csv",header=T)#导入数据>group1<-read.csv("data51a.csv

11、",header=T)#导入数据>group2<-read.csv("data51b.csv",header=T)#导入数据>group1<-group1,-1>group2<-group2,-1> mshapiro.test(t(group1)#RJ用mshapiro.test函数检验数据二元正态性> mshapiro.test(t(group2)wp-value血友病A非携带者0.950.1468血友病A携带者0.970.2535从输出结果可知，两组二元正态性检验的伴随概率分别为0.1468与0.2535,不拒绝

12、原假设(Ho：满足二元正态性)，说明两组数据满足二元正态性。(2)分别计算两组数据的协方差矩阵，是否可以认为两者近似相等?利用cov函数计算协方差矩阵> options(digits=2)> cov(group1)> cov(group2)得到协方差矩阵cov=0.0210.0160.016_0.0240.0150.018、8V-!。.0150.024J利用bartllet检验方差齐性>bartlett.test(data5,-1)Bartletttestofhomogeneityofvariancesdata:datanBartlett'sK-squared=

13、1.4,Df=1p-value=0.2448根据Bartlett检验结果可知，在95%的置信区间内，原假设H0(方差齐性)的伴随概率为0.2448,大于0.05,不拒绝原假设，说明方差齐性，认为两组数据的方差近似相等。(3)对训练样本和新观测合并作散点图，不同的类用不同颜色标识;>library(car)>scatterplot(data5$x1,data5$x2,groups=data5$Group,reg.llne=F,boxplot=F,elllpse=T,levels=0.95,smooth=FALSE,xlab="抗血友病因子活性(x1)”,ylab="

14、;血友病抗原含量(x2)",main="样本双变量分布散点图”,grId=TRUE,col=c("green","red"),legend.title="样本类别”,legend.c00rds="bottomright”,legend.columns=2,cex.main=2)样本双变量分布散点图20_40一-0.6-0.4-0.2抗血友病因子活性(x1)0.0OO20一里含原抗病友血图1散点图分布(4)用lda函数做判别分析，即在协方差矩阵相等的情形下作判别分析；在所给样本信息背景下，将血友病A在样本中的携带概率

15、作为先验信息，进行判别分析，具体如下：>group.lda<-lda(data5,-1,factor(data5$Group)>group.lda组别先验矩阵组均值X1X2血友病A非携带者0.4-0.13-0.078血友病A携带者0.6-0.31-0.006CoefficientsoflineardiscriminantsLD1-98从输出结果先验矩阵(Priorprobabilitiesofgroups)一栏可知，血友病A非携带与携带的先验概率分别为0.4与0.6,即样本中有40%的血友病A的非携带者与60%的血友病A的携带者。由组均值(Groupmean§可知，

16、血友病A非携带者与携带者的Xi与X2均值存在差异，初步判断血友病A非携带者中AHFactivity值略高于携带者，且其AHFantigen水平值也略高于血友病A携带值。同时，用lda方法进行判别分析可得到线性判别函数LD1=-9Xi+8X2.用该线性函数可以计算出每个观察在各组的分类函数值，然后据此分类预测。考虑到血友病A的潜在携带概率可能与样本不一致，查阅相关资料，获得血友病A的潜在携带概率为12499/12500。因此，将其作为先验概率对样本数据重新判别，具体程序和结果如下：>group.ldanew<-lda(data5,-1,factor(data5$Group),prio

17、r=c(12499/12500,1/12500)>group.ldanew组别先验矩阵组均值X1X2血友病A非携带者0.4-0.13-0.078血友病A携带者0.6-0.31-0.006CoefficientsoflineardiscriminantsLD1-98(5)用qda函数做判别分析，即在协方差矩阵不相等的情形下作判别分析;>group.qda<-qda(data5,-1,factor(data5$Group)>group.qda组别先验矩阵组均值X1X2血友病A非携带者0.4-0.13-0.078血友病A携带者0.6-0.31-0.006用qda函数进行判别分

18、析，发现分析结果大致与lda相似，由于qda通过二次判别，并未输出判别函数(6)比较方法(4)和方法(5)的误判率>z2<-predict(group.lda,dim=1)$class>table(z2,data5$Group)表1lda函数判别结果group1212672438>z22<-predict(group.qda,dim=1)$class>table(z22,data5$Group)表2qda函数判别结果group1212652440分别利用lda与qda两个函数对数据进行判别，得到上表。在lda判别方法中,30个血友病A非携带者中有26人被正确

19、判为非携带者，有4人被误判为携带者；45个血友病A携带者中有7人被误判为非携带者，有38人被正确判为携带者，据此可得误判率为14.67%。同样方法可计算qda判别方法的误判率为12%,可见利用qda函数判别分析可降低误判率。实验题目二分析报告：(1)对上表中的数据作散点图，不同的类用不同的颜色标识；>data2<-read.csv("data2.csv",head=T)>colnames(data1)<-c("x1","x2","group")> data21:5,> attac

20、h(data2)> scatterplot(data2$x1,data2$x2,groups=data2$group,smooth=FALSE,reg.line=FALSE,ellipse=TRUE,levels=0.95,diagonal="none",xlab="GPA(x1)”,ylab="GMAT(x2)”,main="样本三变量分布散点图”,legend.title="样本类别”)oR-2.53.03.5样本三变量分布散点图O样本类别。12+3GPA(x1)图2散点图2(2)用Ida函数做判别分析，即在协方差矩阵相等的

21、情形下作判别分析;> attach(data2)> bartlett.test(x1group)> bartlett.test(x2-group)>data2.lda<-lda(data2,-3,factor(data2$group)> data2.ldaBartletttestX1X1p-value0.60.2Bartletttest可以得到两者的p值均大于0.05,所以我们接受原假设，不同组的方差相等分类先验概率组均值X1X2接受0.363.4561不接受0.332.5447彳田0.313.0446从输出结果Priorprobabilitiesofgro

22、ups一栏可知，样本中接受、不接受与待定的概率分别为0.36、0.33与0.31。由Groupmeans可知，接受、不接受与待定的大学期间平均成绩GPA(Xi)和研究生管理能力考试GMAT的成绩(X2)均值存在差异，其中接受待定不接受，说明可以进行判别分析，并符合了分数较高者更容易得到考官青睐这一实际。用lda方法进行判别分析，可得到两个线性判别函数LD1=-5.0088Xi0.0086X2与LD2=1.877Xi0.014X2,其中在进行判另分析过程中，前者方差贡献率为96.7%,后者仅为3.3%。再用该线性函数可以计算出每个观察在各组的分类函数值，然后据此分类预测。(3)用qda函数做判别

23、分析，即在协方差矩阵不相等的情形下作判别分析；>data2.qda<-qda(data2,-3,factor(data2$group)>data2.qda分类先验概率组均值X1X2接受0.363.4561不接受0.332.5447彳田0.313.0446(4)比较方法(2)和方法(3)的误判率;>zl<-predict(data2.lda,dim=2)$class>table(zl,group)表3qda函数判别结果group123128012026133224>zq<-predict(data2.qda,dim=2)$class>tabl

24、e(zq,group)表4qda函数判别结果group123130012027031125分别利用lda与qda两个函数对数据进行判别，得到上表。在lda判别方法中，31个接受的学生中有3个被误判为待定，其余均正确判断；28个未接受的学生中有2个被误判为待定，26个待定的学生中有1个被判为接受，1个被判为未接受，据此计算误判率为8.24%。同样方法计算qda函数判别分析的误判率为4.70%,可见qda函数的判别分析误判率更低。(5)现有一新申请者的GPA为3.21,GMAT成绩为497。请将该观测在(1)的散点图中标出，并分别用方法(2)和方法(3)将其归类？你认为哪一种方法更合适？>n

25、ew<-data.frame(x1=3.21,x2=497)>predict(data2.lda,new)类别接受不接受彳田后验概率0.520.000360.48>predict(data2.qda,new)类别接受不接受彳田后验概率0.920.000450.077从输出结果可知，Ida函数与qda函数均将待估观测判为接受，但后验概率中，一次判别接受和待定的概率相近，用二次判别接受的概率很高，不接受的概率近于0,用二次判别的预测结果更好。(6)观察(1)的散点图中第三类的观测点有无异常值？若有，将该异常值剔除后再对新申请者判别归类，结果有无变化?#把数据化为矩阵形式#化为向量

26、形式#得到异常值的位置#画出盒形图>par(mfrow=c(1,2)>a1<-as.matrix(data260:85,1)> b1<-as.vector(a1)>iny1<-which(b1%in%boxplot.stats(b1)$out)> boxplot(b1)>text(1.1,boxplot.stats(b1)$out,label=paste(异常版”,iny1),col="darkgreen")#观察出奇异点> a2<-as.matrix(data260:85,2)#把数据化为矩阵形式>

27、b2<-as.vector(a2)#化为向量形式> iny2<-which(b2%in%boxplot.stats(b2)$out)>#使用boxplot.stat命令得到异常值的位置> boxplot(b2)#画出盒形图> text(1.1,boxplot.stats(b2)$out,label=paste(异常版”,iny2),col="darkgreen")#观察出奇异点我们通过对画箱型图的方式来找到异常点，并标示在图中。>scatterplot(data2$x1,data2$x2,groups=data2$group,smo

28、oth=FALSE,reg.line=FALSE,ellipse=TRUE,levels=0.95,diagonal="none",xlab="GPA(x1)”,ylab="GMAT(x2)",main="样本三变量分布散点图”,legend.title="样本类别")>points(3.15,313,col="blue",cex=2,pch=1)>points(3.5,402,col="blue",cex=2)样本三变量分布散点翱样本筠IO1A2+3OS89O*

29、重新画散点图，并把根据箱型图方法找到的异常点圈起来，可以发现，通过箱型图找到的异常点与肉眼找到的异常点一致，说明这两点确实为异常点。将两点删除，重新预测。>predict(data3.qda,new)类别接受不接受彳田后验概率0.990.000490.0069用qla预测发现，1的后验概率更大了，说明没有改变预测结果，但加大了预测概率。R程序:实验1data5<-read.csv("data5.csv",header=T)#导入数据data5b<-read.csv("data5b.csv",header=T)#导入数据group1<

30、;-read.csv("data51a.csv",header=T)#导入数据group2<-read.csv("data51b.csv",header=T)#导入数据group1<-group1,-1group2<-group2,-1mshapiro.test(t(group1)#RJ用mshapiro.test函数检验数据二元正态性mshapiro.test(t(group2)options(digits=2)cov(group1)cov(group2)bartlett.test(data5,-1)library(car)scatt

31、erplot(data5$x1,data5$x2,groups=data5$Group,reg.line=F,boxplot=F,ellipse=T,lvels=0.95,smooth=FALSE,xlab="抗血友病因子活性(x1)",ylab="血友病抗原含量(x2)",main="样本双变量分布散点图”,grid=TRUE,col=c("green","red"),legend.title="”,legend.c00rds="bottomright”,legend.columns

32、=2,cex.main=2)group.lda<-lda(data5,-1,factor(data5$Group)group.ldagroup.ldanew<-lda(data5,-1,factor(data5$Group),prior=c(12499/12500,1/12500)group.ldanewz2<-predict(group.lda,dim=1)$classtable(z2,data5$Group)z22<-predict(group.qda,dim=1)$classtable(z22,data5$Group)实验二data2<-read.csv(&

33、quot;data2.csv",head=T)a<-data21:31,-3b<-data232:59,-3c<-data260:85,-3colnames(data2)<-c("x1","x2","group")attach(data2)scatterplot(data2$x1,data2$x2,groups=data2$group,smooth=FALSE,reg.line=FALSEellipse=TRUE,levels=0.95,diagonal="none",xlab="GPA(x1)”,ylab="GMAT(x2)”,main="样本三变量分布散点图”,legend.title="样本类别”)attach(data2)bartlett.test(x1group)bartlett.test(x2group)data2.lda<-lda(data2,-3,factor(data2$group)data2.ldazl<-predict(data2.lda,dim=2)$cla