R语言实验十

实验十判别分析、聚类分析与主成分分析【实验类型】验证性【实验学时】2学时【实验目的】1、掌握判别分析的基本原理及三种重要的判别法的使用条件；2、掌握聚类分析的基本原理及聚类方法的划分方法；3、掌握主成分分析的使用条件及其与逐步回归、岭估计的区别与联系。【实验内容】1、判别分析的计算及其结果解释；2、聚类分析的计算及其结果解释；3、主成分分析的计算及其结果解释。【实验结果】第一部分、课件例题：（第九章1-3节中的所有例题）#例9.1的求解：#按矩阵和因子形式输入数据train-matrix(c(24.8,24.1,26.6,23.5,25.5,27.4,22.1,21.6,22.0,22.8,22.7,21.5,22.1,21.4,-2.0,-2.4,-3.0,-1.9,-2.1,-3.1,-0.7,-1.4,-0.8,-1.6,-1.5,-1.0,-1.2,-1.3),ncol=2)sp-factor(rep(1:2,c(6,8)),labels=c(Have,No))library(MASS);lda.sol-lda(train,sp)#判别分析tst-c(23.5,-1.6);predict(lda.sol,tst)$class#预测table(sp,predict(lda.sol)$class)#回代预测的结果是无春旱；回代的结果说明，原本有6个春旱的年份只判断对了5个。##使用二次判别函数qda.sol-qda(train,sp)predict(qda.sol,tst)$class#预测table(sp,predict(qda.sol)$class)#回代二次判别预测的结果是有春旱，这与前面线性判别的预测结果不同；从回代结果来看，可能是有春早更合理一些，因为二次判别的回代正确率是100%。#例9.2的求解：iris#R自带数据集train-sample(1:150,100)#随机选取100个作为训练样本library(MASS);z1-lda(Species~.,iris,prior=c(1,1,1)/3,subset=train)#Species~.表示对所有自变量；先验概率各为1/3(class1-predict(z1,iris[-train,])$class)#iris[-train,]表示在预测函数中使用其余的50个样本作为待测样本sum(class1==iris$Species[-train])#检验预测结果的准确性结果表明：此次预测结果全部正确（随着每次抽取样本的不同，线性判别预测结果的准确率会有相应的变化）#使用二次判别函数qda()z2-lda(Species~.,iris,prior=c(1,1,1)/3,subset=train)(class2-predict(z2,iris[-train,])$class)#加()表示显示结果sum(class2==iris$Species[-train])#检验预测结果的准确性结果表明：二次判别函数的预测结果全部正确。#例9.3#输入数据,生成距离结构x-c(1,2,6,8,11);dim(x)-c(5,1);d-dist(x)#生成系统聚类hc1-hclust(d,single);hc2-hclust(d,complete)hc3-hclust(d,median);hc4-hclust(d,average)#绘出所有树形结构图,并以2×2的形式绘在一张图上opar-par(mfrow=c(2,2),omi=rep(0,4))plot(hc1,hang=-1);plot(hc2,hang=-1)plot(hc3,hang=-1);plot(hc4,hang=-1)par(opar)结果表明：不同的聚类方法所得的结果基本一样。#例9.4的求解：x-c(1.000,0.846,0.805,0.859,0.473,0.398,0.301,0.382,0.846,1.000,0.881,0.826,0.376,0.326,0.277,0.277,0.805,0.881,1.000,0.801,0.380,0.319,0.237,0.345,0.859,0.826,0.801,1.000,0.436,0.329,0.327,0.365,0.473,0.376,0.380,0.436,1.000,0.762,0.730,0.629,0.398,0.326,0.319,0.329,0.762,1.000,0.583,0.577,0.301,0.277,0.237,0.327,0.730,0.583,1.000,0.539,0.382,0.415,0.345,0.365,0.629,0.577,0.539,1.000)#输入相关矩阵names-c(身高,手臂长,上肢长,下肢长,体重,颈围,胸围,胸宽)r-matrix(x,nrow=8,dimnames=list(names,names))#作系统聚类分析，as.dist()作用是将普通矩阵转化为聚类分析用的距离结构d-as.dist(1-r);hc-hclust(d);dend-as.dendrogram(hc)#写一段小程序,其目的是在绘图命令中调用它,使谱系图更好看nP-list(col=3:2,cex=c(2.0,0.75),pch=21:22,bg=c(lightblue,pink),lab.cex=1.0,lab.col=tomato)addE-function(n){if(!is.leaf(n)){attr(n,edgePar)-list(p.col=plum)attr(n,edgetext)-paste(attr(n,members),members)}}#画出谱系图.de-dendrapply(dend,addE)#对谱系图的所有节点应用函数par(mai=c(0.5,0.5,0.0,0.0))plot(de,nodePar=nP)#绘制谱系图#例9.5的求解：#读取数据X-read.table(F:/文档/大学课程/R语言/ch09/consume.dat)#生成距离结构,作系统聚类d-dist(scale(X))#为消除数据在数量级的影响，对数据作标准化变换method=c(complete,average,centroid,ward.D)for(minmethod){hc-hclust(d,m);class-cutree(hc,k=5)#分成5类print(m);print(sort(class))}#使用sort()排序#使用plot()绘制谱系图，用rect.hclust()将地区分成5类类for(minmethod){hc-hclust(d,m)windows()#建立新窗口显示图形plot(hc,hang=-1)re-rect.hclust(hc,k=5,border=red)print(m);print(re)}结果表明：不同的聚类方法所得的结果有所不同。结果表明：cutree()与rect.hclust()得到的聚类结果是相同的。#例9.6X-read.table(F:/文档/大学课程/R语言/ch09/consume.dat)km-kmeans(scale(X),5,nstart=20)kmsort(km$cluster)#对分类情况排序说明：因为kmeans是随机的确定k个点作为初始点，所以动态聚类每一次的运行结果可能是不相同的。#例9.7的求解：student-read.table(F:/文档/大学课程/R语言/ch09/student_body.dat)#读数据student.pr-princomp(student,cor=TRUE)#用相关矩阵作主成分分析，相当于对数据作标准化变换summary(student.pr,loadings=TRUE)说明：Standarddeviation为主成分的标准差；ProportionofVariance为主成分的贡献率；CumulativeProportion为主成分的累积贡献率；Loadings为载荷矩阵。得到主成分与原始变量的线性关系为：##作预测predict(student.pr)分析：从第1主成分的主成分得分来看：较小得分的样本是25号、3号和5号，说明这几个学生身材魁梧，较大得分的样本是11号、15号和29号，说明这几个学生身材瘦小；从第2主成分的主成分得分来看：较大得分的几个样本是23号、19号和4号，说明这几个学生属于“细高”型，而17号、8号和2号样本的主成分得分较小，说明这几个学生身材属于“矮胖”型。##画出主成分的的碎石图par(mai=c(0.5,0.8,0.4,0.2))screeplot(student.pr,type=lines)#选择折线型结果表明：从第二个主成分之后，图线变化趋于平稳，因此可以选择前两个主成分做分析，与前面的结论一致。##画出双坐标图par(mai=c(.8,.8,.5,.5))biplot(student.pr,scale=0.5)#例9.8的求解：x-scan(F:/文档/大学课程/R语言/ch09/man_body.dat,nmax=120)#sum(1:15)=120,读取120个数据names-scan(F:/文档/大学课程/R语言/ch09/man_body.dat,what=character(0),skip=15)#读取变量名R-matrix(1,nrow=16,ncol=16,dimnames=list(names,names))for(iin2:16){for(jin1:(i-1)){R[i,j]-x[(i-2)*(i-1)/2+j];R[j,i]-R[i,j]}}#生成相关矩阵pr-princomp(covmat=R);load-loadings(pr)#主成份par(mai=c(.8,.8,.2,.1))plot(load[,1:2],pch=19)#画散点图text(load[,1],load[,2],adj=c(-0.4,0.3))#例9.9：对法国经济分析数据(见例7.8)作主成分回归分析。france-data.frame(x1=c(149.3,161.2,171.5,175.5,180.8,190.7,202.1,212.4,226.1,231.9,239.0),x2=c(4.2,4.1,3.1,3.1,1.1,2.2,2.1,5.6,5.0,5.1,0.7),x3=c(108.1,114.8,123.2,126.9,132.1,137.7,146.0,154.1,162.3,164.3,167.6),y=c(15.9,16.4,19.0,19.1,18.8,20.4,22.7,26.5,28.1,27.6,26.3))france.pr-princomp(~x1+x2+x3,data=france,cor=TRUE)summary(france.pr)#作主成份分析结果表明：前两个主成分已达到99%的贡献率，因此选择第1和第2主成分。#例9.9的求解：#主成份回归pre-predict(france.pr)#计算样本的主成分得分france$z1-pre[,1];france$z2-pre[,2]#将第1和第2主成分得分存放在数据框france中lm.sol-lm(y~z1+z2,data=france)#回归分析summary(lm.sol)回归系数和回归方程均通过检验，且效果显著，即得到回归方程：#例9.9的求解：#作变换,得到原坐标下的关系表达式beta-coef(lm.sol);A-loadings(france.pr)[,1:2]x.bar-france.pr$center;x.sd-france.pr$scalecoef-A%*%beta[2:3]/x.sdbeta0-beta[1]-x.bar%*%coefc(be