应用多元统计分析作业

应用多元统计分析作业应用多元统计分析作业应用多元统计分析作业应用多元统计分析作业编制仅供参考审核批准生效日期地址：电话：传真:邮编:多元统计分析实验报告实验课程名称多元统计分析实验项目名称多元统计理论的计算机实现年级2013专业应用统计学学生姓名侯杰成绩理学院实验时间：2015年05月07日学生所在学院：理学院专业：应用统计学班级：01姓名侯杰学号9实验组实验时间指导教师李建军实验项目名称多元统计理论的计算机实现实验目的及要求：目的：熟悉R（或SPSS）软件，掌握多元统计分析中多元正态分布均值向量和协差阵的检验，判别方法，聚类分析，主成分分析，因子分析，相应分析内容。要求：程序要有注释，尽量体现多元统计分析多元正态分布均值向量和协差阵的检验，判别方法，聚类分析，主成分分析，因子分析，相应分析内容内容的基本原理。实验硬件及软件平台：计算机、R、网络实验内容（包括实验具体内容、算法分析、源代码等等）：指导教师意见：签名：年月日代码及运行结果分析均值检验问题重述：某医生观察了16名正常人的24小时动态心电图，分析出早晨3小时各小时的低频心电频谱值（LF）、高频心电频谱值（HF），数据见压缩包，试分析这两个指标的各次重复测定均值向量是否有显著差异。代码如下：<-function(data,alpha={data<("",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<(1-pf(Fstat,3,16+16-3-1));cat("p值=",pvalue,"\n");if(pvalue>#结果输出cat('均值向量不存在差异')elsecat('均值向量存在差异');}运行结果及分析：通过运行程序，我们可以得到如下结果：>()p值=均值向量存在差异即LF与HF这两个指标的各次重复测定均值向量存在显著差异。判别分析问题重述：银行的贷款部门需要判别每个客户的信用好坏（是否未履行还贷责任），以决定是否给予贷款。可以根据贷款申请人的年龄（）、受教育程度（）、现在所从事工作的年数（）、未变更住址的年数（）、收入（）、负债收入比例（）、信用卡债务（）、其它债务（）等来判断其信用情况。数据见压缩包。⑴根据样本资料分别用距离判别法、Bayes判别法和Fisher判别法建立判别函数和判别规则。⑵某客户的如上情况资料为（53，1，9，18，50，，，），对其进行信用好坏的判别。代码如下：#距离判别法<-function(x){data<("",header=T,sep=",");#读取数据G1<-data[1:5,];G2<-data[6: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);#计算判别系数alphaw<-t(alpha)%*%(x-xbar);#构造判别函数if(w>=0)#结果输出cat("该客户属于信用好的一类","\n")elsecat("该客户属于信用坏的一类","\n")}#费希尔判别法<-function(x){data<("",header=T,sep=",");#读取数据G1<-data[1:5,];G2<-data[6: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")elsecat("该客户属于信用坏的一类","\n")}运行结果及分析：注：由于在本题的情形下，距离判别与贝叶斯判别等价，故在此处仅选取距离判别进行编程。距离判别的运行结果：>x<-c(53,1,9,18,50,,,>(x)该客户属于信用好的一类费希尔判别的运行结果：>x<-c(53,1,9,18,50,,,>(x)该客户属于信用好的一类从上面的运行结果可以看出该客户属于信用好的一类，即已履行还贷责任。聚类分析问题重述：下表（数据见压缩包）是某年我国16个地区农民支出情况的抽样调查数据，每个地区调查了反映每人平均生活消费支出情况的六个经济指标。试使用系统聚类法和K均值法对这些地区进行聚类分析，并对结果进行分析比较。代码如下：#系统聚类法data<("",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="");#离差平方和法###绘出四种方法情况下的谱系图和聚类情况opar<-par(mfrow=c(2,2));plot(Clu1,labels=data[,1]);re1<(Clu1,k=5,border="red");box();plot(Clu2,labels=data[,1]);re2<(Clu2,k=5,border="red");box();plot(Clu3,labels=data[,1]);re3<(Clu3,k=5,border="red");box();plot(Clu4,labels=data[,1]);re4<(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],=;axis(2,at=1:5,labels=1:5,=;box();cut2<-cutree(Clu2,k=5);plot(cut2,pch=cut2,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE);axis(1,at=1:16,labels=data[,1],=;axis(2,at=1:5,labels=1:5,=;box();cut3<-cutree(Clu3,k=5);plot(cut3,pch=cut3,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE);axis(1,at=1:16,labels=data[,1],=;axis(2,at=1:5,labels=1:5,=;box();cut4<-cutree(Clu4,k=5);plot(cut4,pch=cut4,ylab="类别编号",xlab="省市",main="聚类的成员",axes=FALSE);axis(1,at=1:16,labels=data[,1],=;axis(2,at=1:5,labels=1:5,=;box();#K均值聚类法data<("",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],=;axis(2,at=1:5,labels=1:5,=;box();legend("topright",c("第一类","第二类","第三类","第四类","第五类"),pch=1:5,cex=;###绘制类中心变量取值折线图plot(Cluk$centers[1,],ylim=c(0,82),xlab="聚类变量",ylab="组均值（类中心）",main="各类聚类变量均值的变化折线图",axes=FALSE);axis(1,at=1:6,labels=c("食品","衣着","燃料","住房","交通和通讯","娱乐教育文化"),=;box();lines(1:6,Cluk$centers[1,],lty=2,col=2);lines(1:6,Cluk$centers[2,],lty=2,col=2);lines(1:6,Cluk$centers[3,],lty=3,col=3);lines(1:6,Cluk$centers[4,],lty=4,col=4);lines(1:6,Cluk$centers[5,],lty=5,col=5);legend("topleft",c("第一类","第二类","第三类","第四类","第五类"),lty=1:5,col=1:5,cex=运行结果及分析：系统聚类法首先输入数据，计算样本间的欧几里得距离，然后用最短距离法，最长距离法，重心法，离差平方和法进行聚类分析，绘出谱系图和直观分类图。分类结果如下：从结果可以看出最短距离法分类的情况如下：第一类：上海第二类：北京第三类：浙江第四类：山东第五类：其余地区最长距离法分类的情况如下：第一类：上海第二类：北京，浙江第三类：吉林，安徽，福建，江西第四类：辽宁，天津，江苏第五类：山西，河北，河南，山东，内蒙，黑龙江重心法分类的情况如下：第一类：上海第二类：北京第三类：浙江第四类：山西，河北，河南第五类：其余地区离差平方和法分类的情况：第一类：上海第二类：北京，浙江第三类：辽宁，天津，江苏第四类：吉林，安徽，福建，江西第五类：其余地区K均值聚类法的结果如下，将结果分成三类分类情况如下：第一类：北京，上海第二类：天津，辽宁，吉林，江苏，浙江，安徽，福建，江西第三类：河北，山西，内蒙，黑龙江，山东，河南并且从第二张图可以看出三类地区的主要差别存在于食品和住房还有交通和通讯的方面，在其他方面则相差无几。主成分分析问题重述：城市用水普及率X1，城市燃气普及率X2，每万人拥有公共交通车辆X3，人均城市道路面积X4，人均公园绿地面积X5，每万人拥有公共厕所X6等六个指标是衡量城市设施水平的主要指标，数据见压缩包，试用主成分分析法对各地区城市设施水平进行综合评价和排序。代码如下：data<("",header=T,sep=",");cormatrix<-cor(data[,2:7]);Result<-eigen(cormatrix);#求相关系数矩阵的特征值和特征向量plot(Result$values,type="b",ylab="特征值（主成分方差）",xlab="特征值编号（主成分编号）");#绘制各主成分的方差变化折线图datamatrix<(data[,2:7]);princomp<-datamatrix%*%Result$vectors;#计算主成分sum<-sum(Result$values);values<(Result$values);contri<-values/sum;#方差贡献率plot(contri,type="b",ylab="方差贡献率",xlab="主成分编号");score<-princomp%*%contri;#计算综合得分place<-data[,1];consequence<(score,place);rank<-c(1:31);final<-consequence[order(-consequence$score),];#对结果进行排序total<-cbind(final,rank);total运行结果及分析：>total>totalscoreplacerank15山东13天津210江苏311浙江413福建52河北619广东714江西89上海91北京1031新疆1130宁夏126辽宁1322重庆1412安徽1520广西1617湖北174山西1827陕西1918湖南2029青海2121海南2223四川237吉林245内蒙古258黑龙江2625云南2728甘肃2816河南2926西藏3024贵州31上面的结果中，rank列为排名，从上到下按得分从大到小排序。此结果可能会有些许疑问，那就是北京上海等强势地区为何排名却少许靠后。从数据方面来看，我们可以看到例如X3每万人拥有公共交通车辆，X4人均城市道路面积，X5人均公园绿地面积，X6每万人拥有公共厕所的数量等这些涉及到人均或者人口数量的指标对于人口十分密集的地区来说得分应该不会太高，因为北京上海这种人口十分密集的地区虽然强势，但是排名却不是十分靠前。因子分析问题重述：利用因子分析方法分析下列30个学生成绩的因子构成，并分析各个学生较适合学文科还是理科。代码如下：library(psych);data<("",header=T,sep=",");correlations<-cor(data);#计算相关系数矩阵fa<-fa(correlations,nfactors=2,rotate="varimax",fm="pa",scores="regression");#设定因子数为2，采用正交旋转进行因子分析###画出因子图形(fa,labels=rownames(fa$loadings));(fa,simple=FALSE);data<(data);faFs<-data%*%fa$weight;#计算因子得分outcome<-c();for(iin1:30){ifelse((faFs[i,1]>faFs[i,2]),outcome<-c(outcome,"文科"),outcome<-c(outcome,"理科"))};outcome<(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""