R-多元统计分析上机讲义

应用多元统计分析

R实验上机讲义TOC\o"1-5"\h\z应用多元统计分析4AppliedMultivariateStatisticalAnalysis4第一章绪论4第二章矩阵4矩阵的建立4矩阵的下标（index）与子集（元素）的提取6矩阵四则运算7矩阵的加减运算7矩阵的相乘8矩阵的求逆8矩阵的其他一些代数运算8求转置矩阵8提取对角元素8矩阵的合并与拉直8方阵的行列式9矩阵的特征根和特征向量9其它函数9矩阵的统计运算11求均值11标准化11减去中位数11第三章多元正态分布及参数的估计12绘制二元正态密度函数及其相应等高线图12多元正态分布的参数估计14多元正态总体的相关量14极大似然估计14第四章多元正态总体参数的假设检验15几个重要统计量的分布15单总体均值向量的检验及置信域16均值向量的检验16样本协方差阵的特征值和特征向量17多总体均值向量的检验17两正态总体均值向量的检验17多个正态总体均值向量的检验-多元方差分析19协方差阵的检验204.4.2多总体协方差阵的检验20独立性检验20正态性检验21第五章判别分析22距离判别22马氏距离22两总体的距离判别22多个总体的距离判别26贝叶斯判别法及广义平方距离判别法26先验用一率（先知知识）26广义平方距离26后验概率（条件概率）27贝叶斯判别准则27费希尔（Fisher）判别29第六章聚类分析30距离和相似系数30距离31数据中心化与标准化变换31相似系数31系统聚类法31类个数的确定34动态聚类法36变量聚类方法36第七章主成分分析37样本的主成分38主成分分析的应用39第八章因子分析42参数估计方法42方差最大的正交旋转45因子得分45第九章对应分析方法46第十章典型相关分析48应用多元统计分析AppliedMultivariateStatisticalAnalysis第一章绪论在实际问题中，很多随机现象涉及到的变量不是一个，而是经常是多个变量，并且这些变量间又存在一定的联系。我们经常需要处理多个变量的观测数据，如果用一元统计方法，由于忽视了各个变量之间可能存在的相关性，一般说来，丢失信息太多，分析的结果不能客观全面反映数据所包含的内容，因此，我们就需要用到多元统计的方法。多元统计分析(MultivariateStatisticalAnalysis)也称多变量统计分析、多因素统计分析或多元分析，是研究客观事物中多变量(多因素或多指标)之间的相互关系和多样品对象之间差异以及以多个变量为代表的多元随机变量之间的依赖和差异的现代统计分析理论和方法。多元统计分析是解决实际问题的有效的数据处理方法。随着电子计算机使用的日益普及，多元统计统计方法已广泛地应用于自然科学、社会科学的各个方面。第二章矩阵矩阵即是二维的数组，它非常的重要，以至于需要单独讨论。由于矩阵应用非常广泛，因此对它定义了一些特殊的应用和操作，R包括许多只对矩阵操作的操作符和函数。2.1矩阵的建立在R中最为常用的是用命令matrix()建立矩阵，而对角矩阵常用函数diag()建立。例如>X<-matrix(1,nr=2,nc=2)X[,1][,2]TOC\o"1-5"\h\z[1,]11[2,]11X<-diag(3)#生成单位阵X[,1][,2][,3][1,]100[2,]010[3,]001diag(2.5,nr=3,nc=5)[,1][,2][,3][,4][,5][1,]2.50.00.000[2,]0.02.50.000[3,]0.00.02.500X<-matrix(1:4,2)#等价于X<-matrix(1:4,2,2)X[,1][,2][1,]13[2,]24rownames(X)<-c("a","b")colnames(X)<-c("c","d")Xcd1324dim(X)22>dimnames(X)[[1]][1]"a""b"[[2]][1]"c""d"注意：①循环准则仍然适用于matrix()，但要求数据项的个数等于矩阵的列数的倍数，否则会出现警告。②矩阵的维数使用c()会得到不同的结果(除非是方阵),因此需要小心。③数据项填充矩阵的方向可通过参数byrow来指定，其缺省是按列填充的(byrow=FALSE),byrow=TRUE|l示按行填充数据。再看几个例子：X<-matrix(1:4,2,4)#按列填充X[,1][,2][,3][,4][1,]1313[2,]2424X<-matrix(1:4,2,3)Warningmessage:Inmatrix(1:4,2,3):数据长度[4]不是矩B^列数[3]的整倍数X<-matrix(1:4,c(2,3))#不经常使用X[,1][,2][1,]13[2,]24X<-matrix(1:4,2,4,byrow=TRUE)#按行填充X[,1][,2][,3][,4][1,]1234[2,]1234因为矩阵是数组的特例，R中数组由函数array()建立，因此矩阵也可以用函数array()来建立，其一般格式为：array(data,dim,dimnames)其中data为一向量，其元素用于构建数组；dim为数组的维数向量(为数值型向量)；dimnames为由各维的名称构成的向量(为字符型向量),缺省为空。看几个例子：A<-array(1:6,c(2,3))A[，1][，2][,3][1,]135[2,]246A<-array(1:4,c(2,3))A[,1][,2][,3][1,]131[2,]242A<-array(1:8,c(2,3))A[,1][,2][,3][1,]135[2,]2462.2矩阵的下标(index)与子集(元素)的提取矩阵的下标可以使用正整数、负整数和逻辑表达式，从而实现子集的提取或修改。考查矩阵x<-matrix(1:6,2,3)x[,1][,2][,3][1,]135[2,]246?提取一个元素x[2,2][1]4?提取若一个或若干个行或列>x[2,2][1]4>x[2,][1]246>x[，2]34x[,2,drop=FALSE][,1][1，]3[2,]4x[,c(2,3),drop=FALSE][，1][，2][1,]35[2,]46?去掉某一个或若干个行与列x[-1,][1]246x[,-2][,1][,2][1,]15[2,]26?添加与替换元素x[,3]<-NAx[,1][,2][,3][1,]13NA[2,]24NAx[is.na(x)]<-1#缺失值用1代替x[,1][,2][,3][1,]131[2,]241矩阵四则运算矩阵也可以进行四则运算(“+”、”—”、“*”、“/"，"”)，分别解释为矩阵对应元素的四则运算。在实际应用中，比较有实际应用的是矩阵的相加，相减，相乘和矩阵的求逆。矩阵的加减运算一般要求矩阵形状完全相同(dim属性完全相同)，矩阵的相乘一般要求一矩阵的列维数与另一矩阵的行维数相同，而矩阵要求逆的话，一般要求它为一方阵。矩阵的加减运算若A,B为两个形状相同的矩阵，两矩阵的和为C,R中表达式为：C<-A+B两矩阵的差为D,R中表达式为：D<-A-B矩阵也可以与数进行加减，A+5表示A中的每个元素加上5。矩阵的相乘操作符％*%用于矩阵相乘。若矩阵A的列数等于矩阵B的行数，矩阵A乘以矩阵B表示为：A%*%B注：X*Y表示两个矩阵的逐元相乘，而不是X和Y的乘积。矩阵的求逆若矩阵A为一方阵，矩阵的逆可以用下面的命令计算：solve(A)o操作符solve()可以用来求解线性方程组：Ax=b,解为solve(A,b)在数学上，用直接求逆的办法解x<-solve(A)%*%b相比solve(A,b)不仅低效而且还有一种潜在的不稳定性。矩阵的其他一些代数运算求转置矩阵转置函数为t(),矩阵X的转置为t(X)。提取对角元素提取对角元的函数为diag()。例如：X<-matrix(1:4,2,2)diag(X)14事实上，diag()的作用依赖于自变量，diag(vector)返回以自变量(向量)为主对角元素的对角矩阵；diag(matrix)返回由矩阵的主对角元素所组成的向量；diag(k)(k为标量)返回k阶单位阵。矩阵的合并与拉直函数cbind()把几个矩阵横向拼成一个大矩阵，这些矩阵行数应该相同；函数rbind()把几个矩阵列向拼成一个大矩阵，这些矩阵列数应该相同。(如果参与合并的矩阵比其它矩阵行数少或列数少，则循环不足后合并。)例如：m1<-matrix。，nr=2,nc=2)m1[,1][,2][1,]11[2,]11m2<-matrix(2,nr=2,nc=2)m2[,1][,2]TOC\o"1-5"\h\z[1,]22[2,]22rbind(m1,m2)[,1][,2][1,]11[2,]11[3,]22[4,]22cbind(m1,m2)[,1][,2][,3][,4][1,]1122[2,]1122方阵的行列式求方阵的行列式使用det()：X<-matrix(1:4,2)>X[,1][,2][1,]13[2,]24>det(X)[1]-2矩阵的特征根和特征向量函数eigen()用来计算矩阵的特征值和特征向量。这个函数的返回值是一个含有values和vectors两个分量的列表。命令A<-eigen(X)>A$values5.3722813-0.3722813$vectors[,1][,2][1,]-0.5657675-0.9093767[2,]-0.82456480.4159736MatrixfacilitesInthefollowingexamples,AandBarematricesandxandbareavectors.OperatororFunctionDescriptionA*BElement-wisemultiplicationA%*%BMatrixmultiplicationA%o%BOuterproduct.AB'crossprod(A,B)crossprod(A)A'BandA'Arespectively.t(A)Transposediag(x)Createsdiagonalmatrixwithelementsofxintheprincipaldiagonaldiag(A)Returnsavectorcontainingtheelementsoftheprincipaldiagonaldiag(k)Ifkisascalar,thiscreatesakxkidentitymatrix.Gofigure.solve(A,b)Returnsvectorxintheequationb=Ax(i.e.,A-1b)solve(A)InverseofAwhereAisasquarematrix.ginv(A)Moore-PenroseGeneralizedInverseofA.ginv(A)requiresloadingtheMASSpackage.y<-eigen(A)y$valaretheeigenvaluesofAy$vecaretheeigenvectorsofAy<-svd(A)SinglevaluedecompositionofA.y$d=vectorcontainingthesingularvaluesofAy$u=matrixwithcolumnscontaintheleftsingularvectorsofAy$v=matrixwithcolumnscontaintherightsingularvectorsofAR<-chol(A)CholeskifactorizationofA.Returnstheuppertriangularfactor,suchthatR'R=A.y<-qr(A)QRdecompositionofA.y$qrhasanuppertrianglethatcontainsthedecompositionandalowertrianglethatcontainsinformationontheQdecomposition.y$rankistherankofA.y$qrauxavectorwhichcontainsadditionalinformationonQ.y$pivotcontainsinformationonthepivotingstrategyused.cbind(A,B,...)Combinematrices(vectors)horizontally.Returnsamatrix.rbind(A,B,...)Combinematrices(vectors)vertically.Returnsamatrix.rowMeans(A)Returnsvectorofrowmeans.rowSums(A)Returnsvectorofrowsums.colMeans(A)Returnsvectorofcolumnmeans.colSums(A)Returnsvectorofcoumnsums..其它函数交叉乘积(crossproduct),函数为crossprod(),crossprod(X,Y)表示一般的内积X'Y,即X的每一列与Y的每一列的内积组成的矩阵；QR分解，函数为qr(),矩阵X的QR分解为X=QRQ为正交阵，R为上三角阵；等等。2.5矩阵的统计运算函数cov()和cor()分别用于计算矩阵的协方差阵和相关系数阵。矩阵的排列是有方向性的，在R中规定矩阵是按列排的，若没有特别说明，函数max(),min(),median(),var(),sd(),sum(),cumsum(),cumprod(),cummax(),cummin()的使用对于矩阵也是按列计算的，但也可以通过选项MARGIN来改变。下面我们要用到对一个对象施加某种运算的函数apply(),其格式为apply(X,MARGIN,FUN)其中X为参与运算的矩阵，FUN为上面的一个函数或“+”、"-”、“*”、“\"(必须放在引号中),MARGIN=1表示按列计算，MARGIN=2表示按行计算。我们还用到sweep()函数，命令sweep(X,MARGIN,STATS,FUN)表示从矩阵X中按MATGIN计算STATS,并从X中除去(sweepout)。求均值m<-matrix(rnorm(n=12),nrow=3)apply(m,MARGIN=1,FUN=mean)#求各行的均值-0.37738650.38641380.2052353>apply(m,MARGIN=2,FUN=mean)#求各列的均值0.33862020.7320669-0.4624578-0.3225460标准化>scale(m,center=T,scale=T)减去中位数row.med<-apply(m,MARGIN=1,FUN=median)sweep(m,MARGIN=1,STATS=row.med,FUN=-”)第三章多元正态分布及参数的估计绘制二元正态密度函数及其相应等高线图书上例2.2.2,5=(l£=Lp=0时的二元正态密度函数及其等高线图：x<-seq(-3,3,by=0.1)y<-xf<-function(x,y,a=1,b=1,r=0){a1=sqrt(a)b1=sqrt(b)d=1-r*rd1=sqrt(d)*a1*b1z=1/(2*pi*d1)*exp((-x*x/a-y*y/b+2*r*x*y/(a1*b1))/(2*d))}z<-outer(x,y,f)#外积函数persp(x,y,z,xlim=range(x),ylim=range(y),zlim=range(z,na.rm=TRUE),theta=30,nticks=5,ticktype="detailed",sub="(t1=(t2=1,p=0时的二元正态密度函数")#密度函数图contour(x,y,z)#等高线图image(x,y,z)#等高线图，实际数据大小用不同色彩表示所得图形为：

-3ct1=o2=1,p=0时的二元正态密度国数-32-10123-3ct1=o2=1,p=0时的二元正态密度国数-32-10123相应等高线图f,把x的任一个元素与f,把x的任一个元素与y的任意一个元素搭配起来作为f的自变量计算得到新的元素值，当函数缺省时表示乘积情况。对参数进行修改，可以绘制任一二元正态密度函数及其相应的等高线图。多元正态分布的参数估计多元正态总体的相关量设观测数据阵为亚]…正"X=iX?Ahi…Xnp?样本均值向量设工㈤=(&.[…,fJ,1=1,2,…，p,则样本均值向量Xn：Xn=1£?=理圜,由2.5.1可得：>Xn<-apply(x,MARGIN=2,mean)或者>ln<-rep(1,n)>Xn<-(ln%*%x)/nXn即为所求样本均值向量。?样本离差阵(交叉乘积阵)样本离差阵A：A=TX-nXtiXn'。>A<-crossprod(x)-2*Xn%*%t(Xn)或者>m<-diag(1,n)-matrix(1,n,n)/n>A<-t(x)%*%m%*%x“口为所求样本离差阵。?样本协方差阵R中求样本协方差阵的函数为cov()。样本数据阵X勺协方差矩阵S即为：>S<-cov(X)?样本相关阵R中求样本协方差阵的函数为cor()。样本数据阵X勺协方差矩阵R即为：>R<-cor(X)极大似然估计极大似然估计法是建立在极大似然原理基础上的一种统计方法。设总体X,其概率密度函数(连续情况)或分布律(离散情况)为P(剜昉，其中B是未知参数(或未知参数向量)。设X1,X2,…，Xn为取自总体X的样本，则似然函数式阶为：L(8)=L(卬⑸…,8)=FELP俏明求使似然函数达到最大的参数日的值，即极大似然估计值。在单参数场合，在R中可以使用函数optimize()求极大似然估计值。optimize()的调用格式如下：optimize(f=,interval=,lower=min(interval),upper=max(interval),maximum=TRUE,tol=.Machine$double.epsA0.25,…)说明：f是似然函数，interval是参数B的取值范围，lower是8的下界，upper是0的上界，maximum=TRUE是求极大值，否则(maximum=FALSE)表示求函数的极小值，tol是表示求值的精确度，…是对f的附加说明。在多参数场合，在R中用函数0Ptim()或者nlm()来求似然函数的极大值，并求相应的极大值点。optim()的调用格式如下：0Ptim(par,fn,gr=NULL,method=c("Nelder-Mead","BFGS","CG","L-BFGS-B","SANN"),lower=-Inf,upper=Inf,control=list(),hessian=FALSE,…)nlm()的定义如下：nlm(f,p,hessian=FALSE,typsize=rep(1,length(p)),fscale=1,print.level=0,ndigit=12,gradtol=1e-6,stepmax=max(1000*sqrt(sum((p/typsize)A2)),1000),steptol=1e-6,iterlim=100,check.analyticals=TRUE,…)三者主要区别是：函数nlm()仅使用牛顿-拉夫逊算法求函数的最小值点；函数0Ptim()提供method选项给出的5种方法中的一种进行优化；上面二个可用于多维函数的极值问题，，而函数optimize()仅适用于一维函数，但可以用于最大与最小值点。(具体选项见帮助。)第四章多元正态总体参数的假设检验在一元统计中，用于检验一元正态总体参数N,小3的抽样分布有好分布，t分布、F分布风，它们都是来自总体NS，口弓的随机样本导出的检验统计量。推广到多元正态总体后，也有相应于以上三个常用分布的统计量：威沙特(Wishart)统计量，霍特林(Hotelling)T?统计量，威尔克斯(Wilks)A统计量，这些统计量是多元统计分析所涉及的假设检验问题的基础。几个重要统计量的分布对于多元正态总体来说，存在几个重要的统计量：威沙特(Wishart)统计量，霍特林（Hotelling）T?统计量，威尔克斯（Wilks）A统计量等，讨论这些统计量的分布是多元统计分析所涉及的假设检验问题的基础。单总体均值向量的检验及置信域均值向量的检验书上例3.2.1,R程序如下>x<-matrix(c(3.7,48.5,9.3,5.7,65.1,8.0,3.8,47.2,10.9,3.2,53.2,12.0,3.1,55.5,9.7,4.6,36.1,7.9,2.4,24.8,14.0,7.2,33.176,6.7,47.4,8.5,5.4,54.1,11.3,3.9,36.9,12.7,4.5,58.8,12.3,3.5,27.8,9.8,4.5,40.2,8.4,1.5,13.5,10.1,8.5,56.4,7.1,4.5,71.6,8.2,6.5,52.8,10.9,4.1,44.1,11.2,5.5,40.9,9.4),20,3,byrow=TRUE)n<-20p<-3u0<-c(4,50,10)#所给总体均值>ln<-rep(1,20)x0<-(ln%*%x)/n#样本均值xm<-x0-u0>mm<-diag(1,20)-matrix(1,20,20)/na<-t(x)%*%mm%*%x#样本离差阵ai=solve(a)dd=xm%*%ai%*%t(xm)d2=(n-1)*ddt2=n*d2;f<-(n-p)*t2/((n-1)*p)#检验统计量f[,1][1,]2.904546fa<-qf(0.95,p,n-p)#自由度为(p,n-p)的F分布的0.95分位数fa[1]3.196777b<-1-pf(f,p,n-p)#尾概率值b[,1][1,]0.06492834>beta<-pf(fa,p,n-p,t2)#犯第二类错误的概率（假设总体均值[1=样本均值m）>beta<-pf(fa,p,n-p,t2)#>beta[1]0.3616381取检验水平为a=0.05,取检验水平为a=0.05,F=2.904546u3.196777=Fa,也可得也相容。在这种情况下，可能犯第二类错误，概率为p=F(F<3.2|p=x0)=0.3616(假定总体均值卜=内学％,取1tl=支0)。样本协方差阵的特征值和特征向量书上例3.2.2,R程序为：x<-matrix(c(3.7,48.5,9.3,5.7,65.1,8.0,3.8,47.2,10.9,3.2,53.2,12.0,3.1,55.5,9.7,4.6,36.1,7.9,2.4,24.8,14.0,7.2,33.1,7.6,6.7,47.4,8.5,5.4,54.1,11.3,3.9,36.9,12.7,4.5,58.8,12.3,3.5,27.8,9.8,4.5,40.2,8.4,1.5,13.5,10.1,8.5,56.4,7.1,4.5,71.6,8.2,6.5,52.8,10.9,4.1,44.1,11.2,5.5,40.9,9.4),20,3,byrow=TRUE)s<-cov(x)s[,1][,2][,3][1,]2.87936810.0100-1.809053[2,]10.010000199.7884-5.640000[3,]-1.809053-5.64003.627658a<-eigen(s)a$values200.4624644.5315911.301392$vectors[,1][,2][,3][1,]-0.05084144-0.573703640.81748351[2,]-0.998283520.05302042-0.02487655[3,]0.029071560.817345080.575414524.3多总体均值向量的检验两正态总体均值向量的检验书上例3.3.1,R程序为：n<-10m<-10p<-4x<-matrix(c(65,75,60,75,70,55,60,65,60,55,+35,50,45,40,30,40,45,40,50,55,25,20,35,40,30,+35,30,25,30,35,60,55,65,70,50,65,60,60,70,75),10)ln<-rep(1,n)x0<-(ln%*%x)/nmx<-diag(1,n)-matrix(1,n,n)/na1<-t(x)%*%mx%*%x>y<-matrix(c(55,50,45,50,55,60,65,50,40,45,+55,60,45,50,50,40,55,60,45,50,40,45,35,50,+30,45,45,35,30,45,65,70,75,70,75,60,75,80,65,70),10)y0<-(ln%*%y)/nmy<-diag(1,n)-matrix(1,n,n)/na2<-t(y)%*%my%*%ya<-a1+a2xy<-x0-y0a