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

共享知识 分享快乐应用多元统计分析R实验上机讲义卑微如蝼蚁、坚强似大象共享知识 分享快乐应用多元统计分析...........................................................................................................................4AppliedMultivariateStatisticalAnalysis...................................................................................4第一章绪论..............................................................................................................................4第二章矩阵..............................................................................................................................42.1矩阵的建立........................................................................................................................42.2矩阵的下标(index)与子集(元素)的提取...........................................................................62.3矩阵四则运算...................................................................................................................72.3.1矩阵的加减运算.....................................................................................................72.3.2矩阵的相乘............................................................................................................82.3.3矩阵的求逆............................................................................................................82.4矩阵的其他一些代数运算.................................................................................................82.4.1求转置矩阵............................................................................................................82.4.2提取对角元素.........................................................................................................82.4.3矩阵的合并与拉直..................................................................................................82.4.4方阵的行列式..........................................................................................................92.4.5矩阵的特征根和特征向量.....................................................................................92.4.6其它函数................................................................................................................92.5矩阵的统计运算..............................................................................................................112.5.1求均值..................................................................................................................112.5.2标准化..................................................................................................................112.5.3减去中位数..........................................................................................................11第三章多元正态分布及参数的估计.....................................................................................123.1绘制二元正态密度函数及其相应等高线图..................................................................123.2多元正态分布的参数估计..............................................................................................143.2.1多元正态总体的相关量.......................................................................................143.2.2极大似然估计.......................................................................................................14第四章多元正态总体参数的假设检验.................................................................................154.1几个重要统计量的分布..................................................................................................154.2单总体均值向量的检验及置信域..................................................................................164.2.1均值向量的检验....................................................................................................164.2.2样本协方差阵的特征值和特征向量....................................................................174.3多总体均值向量的检验...................................................................................................174.3.1两正态总体均值向量的检验...............................................................................174.3.2多个正态总体均值向量的检验-多元方差分析..................................................194.4协方差阵的检验...............................................................................................................204.4.2多总体协方差阵的检验.......................................................................................204.5独立性检验......................................................................................................................204.6正态性检验......................................................................................................................21第五章判别分析.....................................................................................................................225.1距离判别..........................................................................................................................225.1.1马氏距离..............................................................................................................225.1.2两总体的距离判别...............................................................................................225.1.3多个总体的距离判别...........................................................................................265.2贝叶斯判别法及广义平方距离判别法...........................................................................26卑微如蝼蚁、坚强似大象共享知识 分享快乐5.2.1先验概率(先知知识)............................................................................................265.2.2广义平方距离.......................................................................................................265.2.3后验概率（条件概率）.......................................................................................275.2.4贝叶斯判别准则...................................................................................................275.3费希尔（Fisher）判别.....................................................................................................29第六章聚类分析.....................................................................................................................306.2距离和相似系数...............................................................................................................306.2.1距离.......................................................................................................................316.2.2数据中心化与标准化变换....................................................................................316.2.3相似系数...............................................................................................................316.3系统聚类法.....................................................................................................................316.4类个数的确定..................................................................................................................346.5动态聚类法......................................................................................................................366.7变量聚类方法..................................................................................................................36第七章主成分分析.................................................................................................................377.2样本的主成分.................................................................................................................387.3主成分分析的应用..........................................................................................................39第八章因子分析.....................................................................................................................428.3参数估计方法.................................................................................................................428.4方差最大的正交旋转......................................................................................................458.5因子得分.........................................................................................................................45第九章对应分析方法.............................................................................................................46第十章典型相关分析.............................................................................................................48卑微如蝼蚁、坚强似大象共享知识 分享快乐应用多元统计分析AppliedMultivariateStatisticalAnalysis第一章绪论在实际问题中，很多随机现象涉及到的变量不是一个， 而是经常是多个变量， 并且这些变量间又存在一定的联系。我们经常需要处理多个变量的观测数据，如果用一元统计方法，由于忽视了各个变量之间可能存在的相关性， 一般说来，丢失信息太多，分析的结果不能客观全面反映数据所包含的内容，因此，我们就需要用到多元统计的方法。多元统计分析 (MultivariateStatisticalAnalysis) 也称多变量统计分析、多因素统计分析或多元分析，是研究客观事物中多变量 (多因素或多指标 )之间的相互关系和多样品对象之间差异以及以多个变量为代表的多元随机变量之间的依赖和差异的现代统计分析理论和方法。 多元统计分析是解决实际问题的有效的数据处理方法。 随着电子计算机使用的日益普及， 多元统计统计方法已广泛地应用于自然科学、社会科学的各个方面。第二章矩阵矩阵即是二维的数组，它非常的重要，以至于需要单独讨论。由于矩阵应用非常广泛，因此对它定义了一些特殊的应用和操作，R包括许多只对矩阵操作的操作符和函数。2.1矩阵的建立在R中最为常用的是用命令 matrix() 建立矩阵，而对角矩阵常用函数 diag() 建立。例如X<-matrix(1,nr=2,nc=2)X[,1][,2][1,]11[2,]11>X<-diag(3)#生成单位阵>X[,1][,2][,3][1,]100[2,]010[3,]001>diag(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.500卑微如蝼蚁、坚强似大象共享知识 分享快乐X<-matrix(1:4,2)#等价于X<-matrix(1:4,2,2)X[,1][,2][1,] 1 3[2,] 2 4rownames(X)<-c("a","b")colnames(X)<-c("c","d")Xda13b24dim(X)[1]22dimnames(X)[[1]][1]"a""b"[[2]][1]"c""d"注意：①循环准则仍然适用于matrix()，但要求数据项的个数等于矩阵的列数的倍数,否则会出现警告。②矩阵的维数使用 c()会得到不同的结果 (除非是方阵),因此需要小心。③数据项填充矩阵的方向可通过参数 byrow来指定, 其缺省是按列填充的byrow=FALSE)，byrow=TRUE表示按行填充数据。再看几个例子：>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]不是矩阵列数[3]的整倍数>X<-matrix(1:4,c(2,3))#不经常使用>X[,1][,2][1,]13[2,]24>X<-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,] 1 3 1[2,] 2 4 2A<-array(1:8,c(2,3))A[,1][,2][,3][1,] 1 3 5[2,] 2 4 62.2矩阵的下标(index)与子集(元素)的提取矩阵的下标可以使用正整数、 负整数和逻辑表达式， 从而实现子集的提取或修改。 考查矩阵x<-matrix(1:6,2,3)x[,1][,2][,3][1,] 1 3 5[2,] 2 4 6提取一个元素x[2,2][1]4提取若一个或若干个行或列x[2,2]4>x[2,]246卑微如蝼蚁、坚强似大象共享知识 分享快乐x[,2][1]34x[,2,drop=FALSE][,1][1,] 3[2,] 4x[,c(2,3),drop=FALSE][,1][,2][1,] 3 5[2,] 4 6去掉某一个或若干个行与列x[-1,]246>x[,-2][,1][,2][1,]15[2,] 2 6添加与替换元素x[,3]<-NAx[,1][,2][,3][1,]13NA[2,]24NA>x[is.na(x)]<-1#缺失值用1代替>x[,1][,2][,3][1,]131[2,]2412.3矩阵四则运算矩阵也可以进行四则运算(“+”、“-”、“*”、“/”，“^”)，分别解释为矩阵对应元素的四则运算。在实际应用中，比较有实际应用的是矩阵的相加，相减，相乘和矩阵的求逆。矩阵的加减运算一般要求矩阵形状完全相同（dim属性完全相同），矩阵的相乘一般要求一矩阵的列维数与另一矩阵的行维数相同，而矩阵要求逆的话，一般要求它为一方阵。 矩阵的加减运算若A，B为两个形状相同的矩阵，两矩阵的和为 C，R中表达式为：C<-A+B两矩阵的差为 D，R中表达式为：D<-A-B卑微如蝼蚁、坚强似大象共享知识 分享快乐矩阵也可以与数进行加减， A+5表示A中的每个元素加上 5。 矩阵的相乘操作符%*% 用于矩阵相乘。若矩阵 A的列数等于矩阵 B的行数，矩阵 A乘以矩阵表示为：A%*%B注：X*Y表示两个矩阵的逐元相乘，而不是 X和Y的乘积。 矩阵的求逆若矩阵A为一方阵，矩阵的逆可以用下面的命令计算： solve(A)。操作符solve()可以用来求解线性方程组： Ax=b，解为solve(A,b)在数学上，用直接求逆的办法解x<-solve(A)%*%b相比solve(A,b)不仅低效而且还有一种潜在的不稳定性。2.4矩阵的其他一些代数运算 求转置矩阵转置函数为 t() ，矩阵X的转置为 t(X)。 提取对角元素提取对角元的函数为 diag()。例如：X<-matrix(1:4,2,2)diag(X)[1]14事实上，diag()的作用依赖于自变量， diag（vector）返回以自变量（向量）为主对角元素的对角矩阵；diag（matrix）返回由矩阵的主对角元素所组成的向量；diag（k）（k为标量）返回k阶单位阵。矩阵的合并与拉直函数cbind（）把几个矩阵横向拼成一个大矩阵，这些矩阵行数应该相同；函数 rbind（）把几个矩阵列向拼成一个大矩阵，这些矩阵列数应该相同。 （如果参与合并的矩阵比其它矩阵行数少或列数少，则循环不足后合并。）例如：>m1<-matrix(1,nr=2,nc=2)>m1[,1][,2][1,] 1 1[2,] 1 1m2<-matrix(2,nr=2,nc=2)m2卑微如蝼蚁、坚强似大象共享知识 分享快乐[,1][,2][1,] 2 2[2,] 2 2rbind(m1,m2)[,1][,2][1,] 1 1[2,] 1 1[3,] 2 2[4,] 2 2cbind(m1,m2)[,1][,2][,3][,4][1,] 1 1 2 2[2,] 1 1 2 2方阵的行列式求方阵的行列式使用 det()：X<-matrix(1:4,2)>X[,1][,2][1,] 1 3[2,] 2 4det(X)[1]-2 矩阵的特征根和特征向量函数eigen()用来计算矩阵的特征值和特征向量。这个函数的返回值是一个含有values和vectors两个分量的列表。命令A<-eigen(X)>A$values[1]5.3722813-0.3722813$vectors[,1][,2][1,]-0.5657675-0.9093767[2,]-0.82456480.4159736Inthefollowingexamples,AandBarematricesandxandbareavectors.OperatororFunctionDescriptionA*B卑微如蝼蚁、坚强似大象共享知识 分享快乐Element-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=QR，Q为正交阵，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)# 求各行的均值[1]-0.3773865 0.3864138 0.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= -”)卑微如蝼蚁、坚强似大象共享知识 分享快乐第三章多元正态分布及参数的估计3.1绘制二元正态密度函数及其相应等高线图书上例， 时的二元正态密度函数及其等高线图：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=" σ1=σ2=1，ρ=0时的二元正态密度函数 ")密度函数图contour(x,y,z) # 等高线图image(x,y,z) # 等高线图，实际数据大小用不同色彩表示所得图形为：卑微如蝼蚁、坚强似大象共享知识 分享快乐相应等高线图Outer(x,y,f)是一个一般性的外积函数，调用函数f，把x的任一个元素与y的任意一个元素搭配起来作为f的自变量计算得到新的元素值，当函数缺省时表示乘积情况。卑微如蝼蚁、坚强似大象共享知识 分享快乐对参数进行修改，可以绘制任一二元正态密度函数及其相应的等高线图。3.2多元正态分布的参数估计 多元正态总体的相关量设观测数据阵为样本均值向量设 ，=1，2， ，，则样本均值向量 Xn： ,由可得：>Xn<-apply(x,MARGIN=2,mean)或者ln<-rep(1,n)Xn<-(ln%*%x)/nXn即为所求样本均值向量。样本离差阵(交叉乘积阵)样本离差阵A： 。>A<-crossprod(x)-2*Xn%*%t(Xn)或者m<-diag(1,n)-matrix(1,n,n)/nA<-t(x)%*%m%*%xA即为所求样本离差阵。样本协方差阵R中求样本协方差阵的函数为 cov（）。样本数据阵 X的协方差矩阵S即为：>S<-cov(X)样本相关阵R中求样本协方差阵的函数为 cor（）。样本数据阵 X的协方差矩阵R即为:>R<-cor(X) 极大似然估计极大似然估计法是建立在极大似然原理基础上的一种统计方法。设总体 X，其概率密度函数（连续情况）或分布律（离散情况）为 ，其中 是未知参数(或未知参数向量 )。设X1，X2，⋯，Xn为取自总体X的样本，则似然函数 为：卑微如蝼蚁、坚强似大象共享知识 分享快乐⋯, )=求使似然函数达到最大的参数 的值，即极大似然估计值。在单参数场合，在R中可以使用函数optimize()求极大似然估计值。optimize()的调用格式如下：optimize(f=,interval=,lower=min(interval),upper=max(interval),maximum=TRUE,tol= .Machine$double.eps^0.25, ⋯)说明：f是似然函数， interval是参数 的取值范围， lower是 的下界，upper是 的上界，maximum=TRUE是求极大值，否则(maximum=FALSE)表示求函数的极小值，tol是表示求值的精确度，⋯是对f的附加说明。在多参数场合，在R中用函数optim()或者nlm()来求似然函数的极大值，并求相应的极大值点。optim()的调用格式如下：optim(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)^2)),1000),steptol=1e-6,iterlim=100,check.analyticals=TRUE, ⋯)三者主要区别是：函数 nlm()仅使用牛顿-拉夫逊算法求函数的最小值点；函数 optim()提供method选项给出的5种方法中的一种进行优化；上面二个可用于多维函数的极值问题 ,，而函数optimize()仅适用于一维函数，但可以用于最大与最小值点。（具体选项见帮助。）第四章多元正态总体参数的假设检验在一元统计中，用于检验一元正态总体参数 ， 的抽样分布有 分布， 分布、F分布风，它们都是来自总体 的随机样本导出的检验统计量。推广到多元正态总体后，也有相应于以上三个常用分布的统计量：威沙特（ Wishart）统计量，霍特林（ Hotelling ）统计量，威尔克斯（ Wilks） 统计量，这些统计量是多元统计分析所涉及的假设检验问题的基础。4.1几个重要统计量的分布对于多元正态总体来说 ,存在几个重要的统计量 : 威沙特（Wishart）统计量，霍特林卑微如蝼蚁、坚强似大象共享知识 分享快乐（Hotelling ） 统计量，威尔克斯（ Wilks） 统计量等,讨论这些统计量的分布是多元统计分析所涉及的假设检验问题的基础。4.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)>n<-20>p<-3>u0<-c(4,50,10)# 所给总体均值>ln<-rep(1,20)>x0<-(ln%*%x)/n # 样本均值xm<-x0-u0mm<-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.904546>fa<-qf(0.95,p,n-p)# 自由度为(p,n-p) 的F分布的0.95分位数>fa[1]3.196777>b<-1-pf(f,p,n-p)# 尾概率值>b[,1][1,]0.06492834>beta<-pf(fa,p,n-p,t2)# 犯第二类错误的概率（假设总体均值 ）>beta[1]0.3616381取检验水平为 0.05，由尾概率值 p=0.06492834 0.05= ，可得 相容；同样由卑微如蝼蚁、坚强似大象共享知识 分享快乐F=2.904546 3.196777=Fa，也可得 相容。在这种情况下，可能犯第二类错误，概率为=0.3616（假定总体均值 ）。样本协方差阵的特征值和特征向量书上例，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$values[1]200.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多总体均值向量的检验 两正态总体均值向量的检验书上例，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)/n>mx<-diag(1,n)-matrix(1,n,n)/n卑微如蝼蚁、坚强似大象共享知识 分享快乐a1<-t(x)%*%mx%*%xy<-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)/n>my<-diag(1,n)-matrix(1,n,n)/n>a2<-t(y)%*%my%*%y>a<-a1+a2>xy<-x0-y0>ai<-solve(a)>dd<-xy%*%ai%*%t(xy)>d2<-(m+n-2)*dd>t2<-n*m*d2/(n+m)f<-(n+m-1-p)*t2/((n+m-2)*p)pp<-1-pf(f,p,m+n-p-1)x0[,1][,2][,3][,4][1,]644330.563>y0[,1][,2][,3][,4][1,]51.5514070.5>a1[,1][,2][,3][,4][1,]490-170-120.0-245[2,]-17051010.0310[3,]-12010322.5260[4,]-245310260.0510>a2[,1][,2][,3][,4][1,]502.560175-7.5[2,]60.039050195.0[3,]175.050450-100.0[4,]-7.5195-100322.5>d2[,1][1,]5.972499>t2[,1][1,]29.86250>f[,1][1,]6.221353>pp[,1]卑微如蝼蚁、坚强似大象共享知识 分享快乐[1,]0.003705807取检验水平为 0.01，根据尾概率值 p=0.003705807 0.01= ，可得应否定 。 多个正态总体均值向量的检验 -多元方差分析书上例，可利用类似例 或例 的程序进行计算得出结论。下面我们用R自带的manova()函数进行分析。程序如下：x<-read.table("D:/data/d332.txt",header=T)x<-as.matrix(x[,1:4])rate<-factor(gl(3,20),labels=c("group1","group2","group3"))fit<-manova(x~rate)summary.aov(fit)# 对每一个变量进行单因素方差分析summary(fit,test="Wilks")# 使用威尔克斯 统计量程序结果：summary.aov(fit)Responsex1:DfSumSqMeanSqFvalue Pr(>F)rate 239066 19533 8.8780.0004401***Residuals 57125409 2200---Signif.codes:0 ‘***’0.001 ‘**’0.01 ‘*’0.05 ‘.’0.1 ‘’1Responsex2:DfSumSqMeanSqFvaluePr(>F)rate 2 4017 20092.82930.06738.Residuals 5740467 710---Signif.codes:0 ‘***’0.001 ‘**’0.01 ‘*’0.05 ‘.’0.1 ‘’1Responsex3:DfSumSqMeanSqFvaluePr(>F)rate 2 13.43 6.720.18380.8326Residuals 572082.50 36.54Responsex4:DfSumSqMeanSqFvaluePr(>F)rate 2 17.20 8.600.47850.6222Residuals 571024.40 17.97卑微如蝼蚁、坚强似大象共享知识 分享快乐>summary(fit,test="Wilks")Df WilksapproxFnumDfdenDf Pr(>F)rate 20.662123.09069 8 1080.003538**Residuals57---Signif.codes:0 ‘***’0.001 ‘**’0.01 ‘*’0.05 ‘.’0.1 ‘’1结果说明：（1） 取检验水平为 0.01，则对四个指标逐项用一元方差分析方法进行检验，由 p值可得三个组指标间只有第一个指标 有显著差异( =0.0004401)；（2） 取检验水平为 0.01，利用威尔克斯 统计量得到 p=0.003538 0.01，故拒绝原假设，即认为三个组的指标之间有显著差异。4.4协方差阵的检验 多总体协方差阵的检验书上例3.4.1，R程序略（类似例3.2.1或例3.3.1）4.5独立性检验书中例，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)n<-20p<-3x0<-(ln%*%x)/nxm<-x0-u0mm<-diag(1,20)-matrix(1,20,20)/na<-t(x)%*%mm%*%xa0<-det(a)a1<-a[1,1]*a[2,2]*a[3,3]v<-a0/a1b<-n-1.5-(p*p*p-3)/(3*p*p-3*3)df<-0.5*(p*(p+1)-2*3)kc<--b*log(v)p0<-1-pchisq(kc,df)kc9.755514>p0卑微如蝼蚁、坚强似大象共享知识 分享快乐[1]0.02076288取检验水平为 0.05，根据尾概率值 p=0.02076288 0.05= ，可得应否定 ，由R软件所的结果与 SAS软件所的结果一致。4.6正态性检验书中例，R程序为：x<-matrix(c(100,99,96,99,96,75,97,68,76,62,67,34,100,97,100,96,78,97,89,88,84,39,78,37),12)n<-12p<-2ln<-rep(1,n)x0<-(ln%*%x)/ns<-cov(x)si<-solve(s)m<-0for(iin1:n){xx0<-x[i,]-x0dd<-xx0%*%si%*%t(xx0)print(c(i,dd))if(dd<=1.386)m<-m+1}1.0000000.8831922.00000000.77873063.0000000.6965184.0000000.7891365.0000002.1881546.0000002.3848757.00000000.87679298.0000002.0336529.00000000.269104110.0000005.04653111.00000000.789168812.0000005.264147>m7>pp<-m/n>pp0.5833333卑微如蝼蚁、坚强似大象共享知识 分享快乐第五章判别分析判别分析是用于判断样品所属类型的一种统计分析方法。判别分析的目的是对已知归类的数据建立由数值指标构成的归类规则，然后把这样的规则应用到未知归类的样品去归类。在生产、科研和日常生活中经常会遇到如何根据观测到的数据资料对所研究的对象进行判别归类的问题。判别分析问题一般可以如下描述：设有k个维总体，其分布特征已知（如已知分布函数分别为，，或知道来自各个总体的训练样本）。对给定的一个新样品X，判断它来自哪个总体。通常我们先对预先得到的来自这k个总体的若干个样品(称为训练样品)进行检验和归类，来决定相应的判别归类问题是否有意义及误判可能性大小。然后再对给定的一个或几个新的样品，进行判别归类，即决定它(们)自哪个总体。解决这个问题可以有多种途径,下面我们分别讨论几种常用的方法，如距离判别、贝叶斯判别、Fisher判别等。R通用程序：首先我们要用命令>library(MASS)MASS包里的lda()针对线性判别分析。加载MASS宏包，再用函数lda()就可完成判别分析，其基本调用格式如下：lda(formula,data,...,subset,na.action)说明：formula用法为groupsx1+x2+⋯，group表明总体来源，x1，x2，⋯表示分类~指标；subset指明训练样本。具体说明见R帮助。5.1距离判别 马氏距离马氏距离定义：样本X和总体其中 为总体均值向量， 为总体协方差阵。 两总体的距离判别判别准则：其中 = ， 为X到总体的距离。关于两总体距离判别的 R程序(参考薛毅教授的《统计建模与 R软件》一书)：卑微如蝼蚁、坚强似大象共享知识 分享快乐discriminiant.distance<-function(TrnX1,TrnX2,TstX=NULL,var.equal=FALSE){if(is.null(TstX)==TRUE) TstX<-rbind(TrnX1,TrnX2)if(is.vector(TstX)==TRUE) TstX<-t(as.matrix(TstX))elseif(is.matrix(TstX)!=TRUE)TstX<-as.matrix(TstX)if(is.matrix(TrnX1)!=TRUE)TrnX1<-as.matrix(TrnX1)if(is.matrix(TrnX2)!=TRUE)TrnX2<-as.matrix(TrnX2)nx<-nrow(TstX)blong<-matrix(rep(0,nx),nrow=1,byrow=TRUE,dimnames=list("blong",1:nx))mu1<-colMeans(TrnX1);mu2<-colMeans(TrnX2)if(var.equal==TRUE||var.equal==T){S<-var(rbind(TrnX1,TrnX2))w<-mahalanobis(TstX,mu2,S)-mahalanobis(TstX,mu1,S)}else{S1<-var(TrnX1);S2<-var(TrnX2)w<-mahalanobis(TstX,mu2,S2)-mahalanobis(TstX,mu1,S1)}for(iin1:nx){if(w[i]>0)blong[i]<-1elseblong[i]<-2}blong}在程序中，输入变量TrnX1、TrnX2表示训练样本X1，X2，其输入格式是数据框，或矩阵（样本按行输入），输入变量TstX是待测样本，其输入格式是数据框，或矩阵（样本按行输入），或向量（一个待测样本）。如果不输入 TstX（缺省值），则待测样本为两个训练样本之和，即训练样本的回代情况。输入变量var.equal是逻辑变量，var.equal==TRUE表示两个总体的协方差相同；否则（缺省值）为不同。在上述程序中，用到马氏距离函数mahalanobis()，该函数的使用格式为mahalanobis(x,center,cov,inverted=FALSE,...)其中x是样本数据构成的向量或矩阵（ p维），center为样本中心， cov为样本的协方差阵。对于书中例，调用discriminiant.distance()进行判别（假设协方差阵相等）：x<-data.frame(x1=c(13.85,22.31,28.82,15.29,28.79),x2=c(2.79,4.67,4.63,3.54,4.90),x3=c(7.8,12.31,16.18,7.50,16.12),卑微如蝼蚁、坚强似大象共享知识 分享快乐x4=c(49.6,47.8,62.15,43.20,58.10))y<-data.frame(x1=c(2.18,3.85,11.40,3.66,12.10),x2=c(1.06,0.80,0.00,2.42,0.00),x3=c(1.22,4.06,3.50,2.14,5.68),x4=c(20.60,47.10,0.00,15.10,0.00))testx<-rbind(c(8.85,3.38,5.17,26.10),c(28.60,2.40,1.20,127.0),c(20.70,6.70,7.60,30.20),c(7.90,2.40,4.30,33.20),c(3.19,3.20,1.43,9.90),c(12.40,5.10,4.43,24.60),c(16.80,3.40,2.31,31.30),c(15.00,2.70,5.02,64.00))discriminiant.distance(x,y,var.equal=TRUE)blong1111122222discriminiant.distance(x,y,testx,var.equal=TRUE)12345678blong21122111由程序结果可得待判样品 2，3，6，7，8属于含钾盐泉(A盆地)，其余三个属于不含钾盐泉(B盆地)。利用R自带函数lda()，该例的R程序如下：w<-read.table(file="D:/data/disc511.txt",header=T)attach(w)names(w)library(MASS)z<-lda(group~x1+x2+x3+x4)z>pred<-predict(z)$class#predict() 是R内置函数，可以将 lda() 的输出应用于训练样品数据进行预测，从而进行对比。table(pred,group)newdata<-rbind(c(8.85,3.38,5.17,26.10),c(28.60,2.40,1.20,127.0),c(20.70,6.70,7.60,30.20),c(7.90,2.40,4.30,33.20),c(3.19,3.20,1.43,9.90),c(12.40,5.10,4.43,24.60),c(16.80,3.40,2.31,31.30),c(15.00,2.70,5.02,64.00))>newdata<-data.frame(newdata)>predict(z,newdata=newdata)>detach(w)R程序结果:Call:卑微如蝼蚁、坚强似大象共享知识 分享快乐lda(group~x1+x2+x3+x4)Priorprobabilitiesofgroups:B0.50.5Groupmeans:x1 x2 x3 x4A21.8124.10611.98252.17B6.6380.8563.32016.56Coefficientsoflineardiscriminants:LD1x1-0.7794490x2-0.6888651x31.4115135x4-0.1192217grouppredABA50B05$classBAABBAAALevels:AB$posteriorA B11.639701e-039.983603e-0121.000000e+001.932625e-8331.000000e+001.269619e-2048.302424e-029.169758e-0151.190922e-069.999988e-0161.000000e+001.129611e-1071.000000e+001.161894e-2681.000000e+007.135903e-22$xLD11.0536512-31.2985593-7.52868290.3947245卑微如蝼蚁、坚强似大象共享知识 分享快乐2.2416596-3.7639282-9.8136273-8.0017623结果说明：Groupmeans:包含了每组的平均向量；Coefficientsoflineardiscriminants:线性判别系数；列联表表明将训练样品数据代入线性判别函数后的判别结果，两组都没有错判；由$class可以看出8个待判样品，待判样品2，3，6，7，8属于含钾盐泉(A盆地)，其余三个属于不含钾盐泉(B盆地)（与上一程序结果一致）；$posterior给出了后验概率值（具体概念见5.2节）；6)$x给出了线性判别函数的数值。 多个总体的距离判别类似与两个总体的情况，多个总体的情况，按照距离最近的原则对 X进行判别归类时，首先计算样品到各类的马氏 (Mahalanobis)距离，然后进行比较，把待判样品判归距离最小的那个总体。（自编关于多个总体距离判别的 R函数可参考《统计建模与 R软件一书》）。5.2贝叶斯判别法及广义平方距离判别法 先验概率(先知知识)设有k个总体 ，假设事先对所研究的问题有一定的认识，这种认识常用先验概率来描述， 即已知这 k个总体各自出现的概率 （验前概率）为 （显然 ， =1），这组验前概率 称为先验概率。 广义平方距离样品X到总体 (=1, ⋯,k)的广义平方距离 为：,其中是样品X到总体 的马氏距离；其中 为第类的组内样本协方差阵。卑微如蝼蚁、坚强似大象共享知识 分享快乐 后验概率（条件概率）当样品X已知时，它属于 的概率就称为后验概率，一般记为 (或 )。 贝叶斯判别准则几个概念：1.错判概率和错判损失；2.关于先验概率的平均损失。定义5.2.1：设有k个总体：，相应的先验概率为（，=1）。如果有判别法 ，使得 带来的平均损失 达最小，则称判别法 符合贝叶斯判别准则，或称 为贝叶斯判别的解。出于例题需要，学习多总体的 Bayes判别程序（两总体情况参考《统计建模与 R软件一书》）：distinguish.bayes<-function(TrnX,TrnG,p=rep(1,length(levels(TrnG))),TstX=NULL,var.equal=FALSE){if(is.factor(TrnG)==FALSE){mx<-nrow(TrnX);mg<-nrow(TrnG)TrnX<-rbind(TrnX,TrnG)TrnG<-factor(rep(1:2,c(mx,mg)))}if(is.null(TstX)==TRUE)TstX<-TrnXif(is.vector(TstX)==TRUE)TstX<-t(as.matrix(TstX))elseif(is.matrix(TstX)!=TRUE) TstX<-as.matrix(TstX)if(is.matrix(TrnX)!=TRUE) TrnX<-