应用多元统计分析作业(共15页).docx

精选优质文档-倾情为你奉上 多元统计分析 实 验 报 告实验课程名称 多元统计分析 实验项目名称 多元统计理论的计算机实现 年 级 2013 专 业 应用统计学 学生姓名 侯杰 成 绩 理 学 院实验时间： 2015 年 05 月07 日学生所在学院：理学院 专业：应用统计学 班级：姓 名侯杰 学 号4实验组 实验时间 16.05.07指导教师李建军实验项目名称多元统计理论的计算机实现实验目的及要求： 目的：熟悉R（或SPSS）软件，掌握多元统计分析中多元正态分布均值向量和协差阵的检验，判别方法，聚类分析，主成分分析，因子分析，相应分析内容。 要求：程序要有注释，尽量体现多元统计分析多元正态分布均值向量和协差阵的检验，判别方法，聚类分析，主成分分析，因子分析，相应分析内容内容的基本原理。实验硬件及软件平台：计算机、 R、 网络实验内容（包括实验具体内容、算法分析、源代码等等）： 指导教师意见：签名： 年 月 日代码及运行结果分析1、 均值检验问题重述：某医生观察了16名正常人的24小时动态心电图，分析出早晨3小时各小时的低频心电频谱值（LF）、高频心电频谱值（HF），数据见压缩包，试分析这两个指标的各次重复测定均值向量是否有显著差异。代码如下：Tsq.test<-function(data,alpha=0.05) data<-as.matrix(read.table("ch37.csv",header=TRUE,sep=",") #读取数据 xdat<-data,2:4; xbar<-apply(xdat,2,mean); #计算LF指标的均值 ydat<-data,5:7; ybar<-apply(ydat,2,mean); #计算HF指标数据 xcov<-cov(xdat); #计算LF样本协差阵 ycov<-cov(ydat); #计算HF样本协差阵 sinv<-solve(xcov+ycov);#求逆矩阵 Tsq<-(16+16-2)*t(sqrt(16*16/(16+16)*(xbar-ybar)%*%sinv%*%sqrt(16*16/(16+16)*(xbar-ybar); #计算T统计量 Fstat<-(16+16-2)-3+1)/(16+16-2)*3)*Tsq; #计算F统计量 pvalue<-as.numeric(1-pf(Fstat,3,16+16-3-1); cat("p值=",pvalue,"n"); if(pvalue>0.05) #结果输出 cat('均值向量不存在差异') else cat('均值向量存在差异'); 运行结果及分析：通过运行程序，我们可以得到如下结果： > Tsq.test() p值= 1.e-14 均值向量存在差异即LF与HF这两个指标的各次重复测定均值向量存在显著差异。2、 判别分析问题重述：银行的贷款部门需要判别每个客户的信用好坏（是否未履行还贷责任），以决定是否给予贷款。可以根据贷款申请人的年龄（）、受教育程度（）、现在所从事工作的年数（）、未变更住址的年数（）、收入（）、负债收入比例（）、信用卡债务（）、其它债务（）等来判断其信用情况。数据见压缩包。根据样本资料分别用距离判别法、Bayes判别法和Fisher判别法建立判别函数和判别规则。某客户的如上情况资料为（53，1，9，18，50，11.20，2.02，3.58），对其进行信用好坏的判别。 代码如下：#距离判别法discrim.dist<-function(x) data<-read.csv("ch49.csv",header=T,sep=","); #读取数据 G1<-data1:5,; G2<-data6:10,; u1<-apply(G1,2,mean); #计算信用好的样本数据均值 u2<-apply(G2,2,mean); #计算信用不好的样本数据均值 s1<-cov(G1); s2<-cov(G2); s<-s1+s2; xbar<-(u1+u2)/2; alpha<-solve(s)%*%(u1-u2); #计算判别系数alpha w<-t(alpha)%*%(x-xbar); #构造判别函数 if(w>=0) #结果输出 cat("该客户属于信用好的一类","n") else cat("该客户属于信用坏的一类","n")#费希尔判别法fisher.test<-function(x) data<-read.csv("ch49.csv",header=T,sep=","); #读取数据 G1<-data1:5,; G2<-data6:10,; n1<-nrow(G1); n2<-nrow(G2); u1<-apply(G1,2,mean); #计算信用好的一组的数据均值 u2<-apply(G2,2,mean); #计算信用不好的一组的样本数据均值 s1<-cov(G1); s2<-cov(G2); E<-s1+s2; B<-n1*n2*(u1-u2)%*%t(u1-u2)/(n1+n1); alpha<-eigen(solve(E)%*%B); vector<-alpha$vectors,1; #提取费希尔判别函数系数 d1<-abs(t(vector)%*%x-t(vector)%*%u1); #计算样本到第一组的费希尔判别函数值 d2<-abs(t(vector)%*%x-t(vector)%*%u2); #计算样本到第二组的费希尔判别函数值 if(d1<d2) #结果输出 cat("该客户属于信用好的一类","n") else cat("该客户属于信用坏的一类","n")运行结果及分析：注：由于在本题的情形下，距离判别与贝叶斯判别等价，故在此处仅选取距离判别进行编程。距离判别的运行结果：> x<-c(53,1,9,18,50,11.20,2.02,3.58)> discrim.dist(x)该客户属于信用好的一类费希尔判别的运行结果：> x<-c(53,1,9,18,50,11.20,2.02,3.58)> fisher.test(x)该客户属于信用好的一类 从上面的运行结果可以看出该客户属于信用好的一类，即已履行还贷责任。3、 聚类分析问题重述：下表（数据见压缩包）是某年我国16个地区农民支出情况的抽样调查数据，每个地区调查了反映每人平均生活消费支出情况的六个经济指标。试使用系统聚类法和K均值法对这些地区进行聚类分析，并对结果进行分析比较。代码如下：#系统聚类法data<-read.csv("ch58.csv",header=T,sep=","); #读取数据Cludata<-data,2:7;Dismatrix<-dist(Cludata,method="euclidean"); #计算样本间的欧几里得距离Clu1<-hclust(d=Dismatrix,method="single"); #最短距离法Clu2<-hclust(d=Dismatrix,method="complete"); #最长距离法Clu3<-hclust(d=Dismatrix,method="centroid"); #重心法Clu4<-hclust(d=Dismatrix,method="ward.D"); #离差平方和法#绘出四种方法情况下的谱系图和聚类情况opar<-par(mfrow=c(2,2);plot(Clu1,labels=data,1);re1<-rect.hclust(Clu1,k=5,border="red");box();plot(Clu2,labels=data,1);re2<-rect.hclust(Clu2,k=5,border="red");box();plot(Clu3,labels=data,1);re3<-rect.hclust(Clu3,k=5,border="red");box();plot(Clu4,labels=data,1);re4<-rect.hclust(Clu4,k=5,border="red");box();par(opar);#绘出直观的分类情形opar<-par(mfrow=c(2,2),las=2);cut1<-cutree(Clu1,k=5);plot(cut1,pch=cut1,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE);axis(1,at=1:16,labels=data,1,cex.axis=0.6);axis(2,at=1:5,labels=1:5,cex.axis=0.6);box();cut2<-cutree(Clu2,k=5);plot(cut2,pch=cut2,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE);axis(1,at=1:16,labels=data,1,cex.axis=0.6);axis(2,at=1:5,labels=1:5,cex.axis=0.6);box();cut3<-cutree(Clu3,k=5);plot(cut3,pch=cut3,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE);axis(1,at=1:16,labels=data,1,cex.axis=0.6);axis(2,at=1:5,labels=1:5,cex.axis=0.6);box();cut4<-cutree(Clu4,k=5);plot(cut4,pch=cut4,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE);axis(1,at=1:16,labels=data,1,cex.axis=0.6);axis(2,at=1:5,labels=1:5,cex.axis=0.6);box();#K均值聚类法data<-read.csv("ch58.csv",header=T,sep=","); #读取数据Cludata<-data,2:7;Cluk<-kmeans(x=Cludata,centers=5,nstart=5); #用K均值聚类分成五类par(mfrow=c(2,1),las=2);cluster<-Cluk$cluster; #保存聚类解plot(cluster,pch=cluster,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE); #绘制各省市聚类解得序列图axis(1,at=1:16,labels=data,1,cex.axis=0.6);axis(2,at=1:5,labels=1:5,cex.axis=0.6);box();legend("topright",c("第一类","第二类","第三类","第四类","第五类"),pch=1:5,cex=0.3);#绘制类中心变量取值折线图plot(Cluk$centers1,ylim=c(0,82),xlab="聚类变量",ylab="组均值（类中心）",main="各类聚类变量均值的变化折线图",axes=FALSE);axis(1,at=1:6,labels=c("食品","衣着","燃料","住房","交通和通讯","娱乐教育文化"),cex.axis=0.6);box();lines(1:6,Cluk$centers1,lty=2,col=2);lines(1:6,Cluk$centers2,lty=2,col=2);lines(1:6,Cluk$centers3,lty=3,col=3);lines(1:6,Cluk$centers4,lty=4,col=4);lines(1:6,Cluk$centers5,lty=5,col=5);legend("topleft",c("第一类","第二类","第三类","第四类","第五类"),lty=1:5,col=1:5,cex=0.2)运行结果及分析：系统聚类法首先输入数据，计算样本间的欧几里得距离，然后用最短距离法，最长距离法，重心法，离差平方和法进行聚类分析，绘出谱系图和直观分类图。分类结果如下：从结果可以看出最短距离法分类的情况如下：第一类：上海第二类：北京第三类：浙江第四类：山东第五类：其余地区最长距离法分类的情况如下：第一类：上海第二类：北京，浙江第三类：吉林，安徽，福建，江西第四类：辽宁，天津，江苏第五类：山西，河北，河南，山东，内蒙，黑龙江重心法分类的情况如下：第一类：上海第二类：北京第三类：浙江第四类：山西，河北，河南第五类：其余地区离差平方和法分类的情况：第一类：上海第二类：北京，浙江第三类：辽宁，天津，江苏第四类：吉林，安徽，福建，江西第五类：其余地区K均值聚类法的结果如下，将结果分成三类分类情况如下：第一类：北京，上海第二类：天津，辽宁，吉林，江苏，浙江，安徽，福建，江西第三类：河北，山西，内蒙，黑龙江，山东，河南 并且从第二张图可以看出三类地区的主要差别存在于食品和住房还有交通和通讯的方面，在其他方面则相差无几。 4、 主成分分析问题重述：城市用水普及率X1，城市燃气普及率X2，每万人拥有公共交通车辆X3，人均城市道路面积X4，人均公园绿地面积X5，每万人拥有公共厕所X6等六个指标是衡量城市设施水平的主要指标，数据见压缩包，试用主成分分析法对各地区城市设施水平进行综合评价和排序。代码如下：data<-read.csv("ch68.csv",header=T,sep=",");cormatrix<-cor(data,2:7);Result<-eigen(cormatrix); #求相关系数矩阵的特征值和特征向量plot(Result$values,type="b",ylab="特征值（主成分方差）",xlab="特征值编号（主成分编号）"); #绘制各主成分的方差变化折线图datamatrix<-as.matrix(data,2:7);princomp<-datamatrix%*%Result$vectors; #计算主成分sum<-sum(Result$values); values<-as.vector(Result$values);contri<-values/sum; #方差贡献率plot(contri,type="b",ylab="方差贡献率",xlab="主成分编号");score<-princomp%*%contri; #计算综合得分place<-data,1;consequence<-data.frame(score,place);rank<-c(1:31); final<-consequenceorder(-consequence$score),; #对结果进行排序total<-cbind(final,rank);total运行结果及分析：> total> total score place rank15 48.17886 山东 13 47.02848 天津 210 46.98910 江苏 311 45.69279 浙江 413 45.67545 福建 52 45.49443 河北 619 45.42103 广东 714 44.84118 江西 89 44.47630 上海 91 44.19970 北京 1031 43.89986 新疆 1130 43.89648 宁夏 126 43.85077 辽宁 1322 43.84919 重庆 1412 43.81336 安徽 1520 43.65408 广西 1617 43.41939 湖北 174 42.94670 山西 1827 42.73741 陕西 1918 41.22774 湖南 2029 40.63502 青海 2121 40.16107 海南 2223 40.02236 四川 237 39.94659 吉林 245 39.33251 内蒙古 258 39.01378 黑龙江 2625 38.63829 云南 2728 37.48472 甘肃 2816 36.96833 河南 2926 36.29778 西藏 3024 35.74812 贵州 31上面的结果中，rank列为排名，从上到下按得分从大到小排序。此结果可能会有些许疑问，那就是北京上海等强势地区为何排名却少许靠后。从数据方面来看，我们可以看到例如X3每万人拥有公共交通车辆，X4人均城市道路面积，X5人均公园绿地面积，X6每万人拥有公共厕所的数量等这些涉及到人均或者人口数量的指标对于人口十分密集的地区来说得分应该不会太高，因为北京上海这种人口十分密集的地区虽然强势，但是排名却不是十分靠前。5、 因子分析问题重述：利用因子分析方法分析下列30个学生成绩的因子构成，并分析各个学生较适合学文科还是理科。代码如下： library(psych); data<-read.csv("ch77.csv",header=T,sep=","); correlations<-cor(data); #计算相关系数矩阵 fa<-fa(correlations,nfactors=2,rotate="varimax",fm="pa",scores="regression"); #设定因子数为2，采用正交旋转进行因子分析 #画出因子图形 factor.plot(fa,labels=rownames(fa$loadings); fa.diagram(fa,simple=FALSE); data<-as.matrix(data); faFs<-data%*%fa$weight; #计算因子得分 outcome<-c(); for(i in 1:30) ifelse(faFsi,1>faFsi,2),outcome<-c(outcome,"文科"),outcome<-c(outcome,"理科") ; outcome<-as.vector(outcome); #根据得分决定文理科 result<-cbind(data,outcome);运行结果及分析：从图上我们看出数学物理化学在第二个因子上载荷较大，语文英语历史在第一个因子载荷上较大，因此我们也就可以将第一个因子和第二个因子理解成日常中的文科和理科接下来根据每个学生的因子得分比较来分析该学生适合文科还是理科，结果如下：> result 数学 物理 化学 语文 历史 英语 outcome 1, "65" "61" "72" "84" "81" "79" "文科" 2, "77" "77" "76" "64" "70" "55" "文科" 3, "67" "63" "49" "65" "67" "57" "文科" 4, "80" "69" "75" "74" "74" "63" "文科" 5, "74" "70" "80" "84" "81" "74" "文科" 6, "78" "84" "75" "62" "71" "64" "文科" 7, "66" "71" "67" "52" "65" "57" "文科" 8, "77" "71" "57" "72" "86" "71" "文科" 9, "83" "100" "79" "41" "67" "50" "理科" 10, "86" "94" "97" "51" "63" "55" "理科" 11, "74" "80" "88" "64" "73" "66" "文科" 12, "67" "84" "53" "58" "66" "56" "文科" 13, "81" "62" "69" "56" "66" "52" "理科" 14, "71" "64" "94" "52" "61" "52" "理科" 15, "78" "96" "81" "80" "89" "76" "文科" 16, "69" "56" "67" "75" "94" "80" "文科" 17, "77" "90" "80" "68" "66" "60" "文科" 18, "84" "67" "75" "60" "70" "63" "文科" 19, "62" "67" "83" "71" "85" "77" "文科" 20, "74" "65" "75" "72" "90" "73" "文科" 21, "91" "74" "97" "62" "71" "66" "理科" 22, "72" "87" "72" "79" "83" "76" "文科" 23, "82" "70" "83" "68" "77" "85" "文科" 24, "63" "70" "60" "91" "85" "82" "文科" 25, "74" "79" "95" "59" "74" "59" "理科" 26, "66" "61" "77" "62" "73" "64" "文科" 27, "90" "82" "98" "47" "71" "60" "理科" 28, "77" "90" "85" "68" "73" "76" "文科" 29, "91" "82" "84" "54" "62" "60" "理科" 30, "78" "84" "100" "51" "60" "60" "理科" 6 相应分析 问题重述：费希尔研究头发颜色与眼睛颜色的关系，抽查了5387人的资料如下表，试对其进行相应分析。数据见压缩包。 代码如下：library(ca);data<-read.csv("ch85.csv",header=T);data1<-as.matrix(data,2:6);rownames(data1)<-data,1;ca<-ca(data1);plot(ca);运行结果及分析：> ca Principal inertias (eigenvalues): 1 2 3 Value 0. 0. 0.Percentage 83.34% 12.87% 3.79% Rows: 蓝色 淡蓝 浅蓝 深蓝Mass 0. 0. 0. 0.ChiDist 0. 0. 0. 0.Inertia 0. 0. 0. 0.Dim. 1 -0. -0. 0. 1.Dim. 2 -1. -0. 1. -0.Columns: 金黄色 红色 褐色 深红 黑色Mass 0. 0. 0. 0. 0.ChiDist 0. 0. 0. 0. 1.Inertia 0. 0. 0. 0. 0.Dim. 1 -1. -0. -0. 1. 2.Dim. 2 -1. 0. 1. -0. -1.从结果我们可以看出眼睛颜色与头发颜色这两个变量之间的关系为：淡蓝色和蓝色的眼睛与金黄色的头发有较大关联；浅蓝色的眼睛和红色还有褐色的头发有较大的关联；深蓝色的眼睛与深绿还有黑色的头发有较大的关联。专心-专注-专业