R-多元统计分析上机讲义.doc
多元分析R实验上机讲义应用多元统计分析R实验上机讲义应用多元统计分析4Applied Multivariate Statistical Analysis4第一章绪论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 矩阵的统计运算102.5.1 求均值102.5.2 标准化102.5.3 减去中位数10第三章多元正态分布及参数的估计103.1 绘制二元正态密度函数及其相应等高线图103.2 多元正态分布的参数估计123.2.1 多元正态总体的相关量123.2.2 极大似然估计13第四章多元正态总体参数的假设检验144.1 几个重要统计量的分布144.2 单总体均值向量的检验及置信域144.2.1均值向量的检验144.2.2样本协方差阵的特征值和特征向量154.3多总体均值向量的检验164.3.1 两正态总体均值向量的检验164.3.2 多个正态总体均值向量的检验-多元方差分析174.4协方差阵的检验184.4.2 多总体协方差阵的检验184.5独立性检验184.6正态性检验19第五章判别分析205.1距离判别215.1.1 马氏距离215.1.2 两总体的距离判别215.1.3 多个总体的距离判别245.2贝叶斯判别法及广义平方距离判别法245.2.1 先验概率(先知知识)245.2.2 广义平方距离255.2.3 后验概率（条件概率）255.2.4 贝叶斯判别准则255.3费希尔（Fisher）判别28第六章聚类分析296.2距离和相似系数296.2.1距离296.2.2数据中心化与标准化变换296.2.3相似系数306.3 系统聚类法306.4类个数的确定326.5动态聚类法346.7变量聚类方法34第七章主成分分析357.2 样本的主成分367.3 主成分分析的应用37第八章因子分析408.3 参数估计方法408.4 方差最大的正交旋转438.5 因子得分43第九章对应分析方法44第十章典型相关分析46应用多元统计分析Applied Multivariate Statistical Analysis第一章 绪论在实际问题中，很多随机现象涉及到的变量不是一个，而是经常是多个变量，并且这些变量间又存在一定的联系。我们经常需要处理多个变量的观测数据，如果用一元统计方法，由于忽视了各个变量之间可能存在的相关性，一般说来，丢失信息太多，分析的结果不能客观全面反映数据所包含的内容，因此，我们就需要用到多元统计的方法。多元统计分析(Multivariate Statistical Analysis)也称多变量统计分析、多因素统计分析或多元分析，是研究客观事物中多变量(多因素或多指标)之间的相互关系和多样品对象之间差异以及以多个变量为代表的多元随机变量之间的依赖和差异的现代统计分析理论和方法。多元统计分析是解决实际问题的有效的数据处理方法。随着电子计算机使用的日益普及，多元统计统计方法已广泛地应用于自然科学、社会科学的各个方面。第二章 矩阵 矩阵即是二维的数组，它非常的重要，以至于需要单独讨论。由于矩阵应用非常广泛，因此对它定义了一些特殊的应用和操作，R 包括许多只对矩阵操作的操作符和函数。2.1矩阵的建立在R中最为常用的是用命令matrix( )建立矩阵，而对角矩阵常用函数diag( )建立。例如> X<-matrix(1,nr=2,nc=2)> X ,1 ,21, 1 12, 1 1> X <- diag(3) # 生成单位阵> X ,1 ,2 ,31, 1 0 02, 0 1 03, 0 0 1> diag(2.5, nr = 3, nc = 5) ,1 ,2 ,3 ,4 ,51, 2.5 0.0 0.0 0 02, 0.0 2.5 0.0 0 03, 0.0 0.0 2.5 0 0> X <- matrix(1:4, 2) # 等价于X <- matrix(1:4, 2, 2)> X ,1 ,21, 1 32, 2 4> rownames(X) <- c("a", "b")> colnames(X) <- c("c", "d")> X c da 1 3b 2 4> dim(X)1 2 2> dimnames(X)11 "a" "b"21 "c" "d"注意：循环准则仍然适用于matrix( )，但要求数据项的个数等于矩阵的列数的倍数, 否则会出现警告。矩阵的维数使用c( )会得到不同的结果(除非是方阵), 因此需要小心。数据项填充矩阵的方向可通过参数byrow来指定, 其缺省是按列填充的（byrow=FALSE)， byrow=TRUE表示按行填充数据。再看几个例子：> X <- matrix(1:4, 2, 4) # 按列填充> X ,1 ,2 ,3 ,41, 1 3 1 32, 2 4 2 4> X <- matrix(1:4, 2, 3)Warning message:In matrix(1:4, 2, 3) : 数据长度4不是矩阵列数3的整倍数> X <- matrix(1:4, c(2, 3) # 不经常使用> X ,1 ,21, 1 32, 2 4> X <- matrix(1:4, 2, 4, byrow=TRUE) # 按行填充> X ,1 ,2 ,3 ,41, 1 2 3 42, 1 2 3 4 因为矩阵是数组的特例,R中数组由函数array( )建立, 因此矩阵也可以用函数array( )来建立，其一般格式为：> array(data, dim, dimnames)其中data为一向量，其元素用于构建数组；dim为数组的维数向量（为数值型向量）；dimnames为由各维的名称构成的向量（为字符型向量)， 缺省为空。看几个例子：> A<-array(1:6,c(2,3)> A ,1 ,2 ,31, 1 3 52, 2 4 6> A<-array(1:4,c(2,3)> A ,1 ,2 ,31, 1 3 12, 2 4 2> A<-array(1:8,c(2,3)> A ,1 ,2 ,31, 1 3 52, 2 4 62.2矩阵的下标(index)与子集(元素)的提取矩阵的下标可以使用正整数、负整数和逻辑表达式，从而实现子集的提取或修改。考查矩阵> x <- matrix(1:6, 2, 3)> x ,1 ,2 ,31, 1 3 52, 2 4 6 提取一个元素> x2,21 4 提取若一个或若干个行或列> x2,21 4> x2,1 2 4 6> x,21 3 4> x,2,drop=FALSE ,11, 32, 4> x,c(2,3),drop=FALSE ,1 ,21, 3 52, 4 6 去掉某一个或若干个行与列> x-1,1 2 4 6> x,-2 ,1 ,21, 1 52, 2 6 添加与替换元素> x,3 <- NA> x ,1 ,2 ,31, 1 3 NA2, 2 4 NA> xis.na(x) <- 1 # 缺失值用1代替> x ,1 ,2 ,31, 1 3 12, 2 4 12.3 矩阵四则运算矩阵也可以进行四则运算(“+”、“-”、“*”、“/”，“”)，分别解释为矩阵对应元素的四则运算。在实际应用中，比较有实际应用的是矩阵的相加，相减，相乘和矩阵的求逆。矩阵的加减运算一般要求矩阵形状完全相同（dim属性完全相同），矩阵的相乘一般要求一矩阵的列维数与另一矩阵的行维数相同，而矩阵要求逆的话，一般要求它为一方阵。2.3.1 矩阵的加减运算若A，B为两个形状相同的矩阵，两矩阵的和为C，R中表达式为：C<-A+B两矩阵的差为D，R中表达式为：D<-A-B矩阵也可以与数进行加减，A+5表示A中的每个元素加上5。2.3.2 矩阵的相乘操作符%*% 用于矩阵相乘。若矩阵A的列数等于矩阵B的行数，矩阵A乘以矩阵B表示为：A%*%B注：X*Y表示两个矩阵的逐元相乘，而不是X和Y的乘积。2.3.3 矩阵的求逆若矩阵A为一方阵，矩阵的逆可以用下面的命令计算：solve(A)。操作符solve( )可以用来求解线性方程组：Ax=b，解为solve(A,b)在数学上，用直接求逆的办法解x <- solve(A) %*% b相比solve(A,b)不仅低效而且还有一种潜在的不稳定性。2.4矩阵的其他一些代数运算2.4.1 求转置矩阵转置函数为t( )，矩阵X的转置为t(X)。2.4.2 提取对角元素提取对角元的函数为diag( )。例如：> X <- matrix(1:4, 2, 2)> diag(X)1 1 4事实上，diag( )的作用依赖于自变量，diag（vector）返回以自变量（向量）为主对角元素的对角矩阵；diag（matrix）返回由矩阵的主对角元素所组成的向量；diag（k）（k为标量）返回k阶单位阵。2.4.3矩阵的合并与拉直函数cbind（）把几个矩阵横向拼成一个大矩阵，这些矩阵行数应该相同；函数rbind（）把几个矩阵列向拼成一个大矩阵，这些矩阵列数应该相同。（如果参与合并的矩阵比其它矩阵行数少或列数少，则循环不足后合并。）例如：> m1 <- matrix(1, nr = 2, nc = 2)> m1 ,1 ,21, 1 12, 1 1> m2 <- matrix(2, nr = 2, nc = 2)> m2 ,1 ,21, 2 22, 2 2> rbind(m1, m2) ,1 ,21, 1 12, 1 13, 2 24, 2 2> cbind(m1, m2) ,1 ,2 ,3 ,41, 1 1 2 22, 1 1 2 22.4.4方阵的行列式求方阵的行列式使用det( )：X<-matrix(1:4, 2)> X ,1 ,21, 1 32, 2 4> det(X)1 -22.4.5 矩阵的特征根和特征向量 函数eigen( ) 用来计算矩阵的特征值和特征向量。这个函数的返回值是一个含有values 和vectors 两个分量的列表。命令A<-eigen(X)> A$values1 5.3722813 -0.3722813$vectors ,1 ,21, -0.5657675 -0.90937672, -0.8245648 0.41597362.4.6 Matrix facilites In the following examples, A and B are matrices and x and b are a vectors.Operator or FunctionDescriptionA * BElement-wise multiplicationA %*% BMatrix multiplicationA %o% BOuter product. AB'crossprod(A,B)crossprod(A)A'B and A'A respectively.t(A)Transposediag(x)Creates diagonal matrix with elements of x in the principal diagonaldiag(A)Returns a vector containing the elements of the principal diagonaldiag(k)If k is a scalar, this creates a k x k identity matrix. Go figure.solve(A, b)Returns vector x in the equation b = Ax (i.e., A-1b)solve(A)Inverse of A where A is a square matrix.ginv(A)Moore-Penrose Generalized Inverse of A. ginv(A) requires loading the MASS package.y<-eigen(A)y$val are the eigenvalues of Ay$vec are the eigenvectors of Ay<-svd(A)Single value decomposition of A.y$d = vector containing the singular values of Ay$u = matrix with columns contain the left singular vectors of A y$v = matrix with columns contain the right singular vectors of AR <- chol(A)Choleski factorization of A. Returns the upper triangular factor, such that R'R = A.y <- qr(A)QR decomposition of A. y$qr has an upper triangle that contains the decomposition and a lower triangle that contains information on the Q decomposition.y$rank is the rank of A. y$qraux a vector which contains additional information on Q. y$pivot contains information on the pivoting strategy used.cbind(A,B,.)Combine matrices(vectors) horizontally. Returns a matrix.rbind(A,B,.)Combine matrices(vectors) vertically. Returns a matrix.rowMeans(A)Returns vector of row means.rowSums(A)Returns vector of row sums.colMeans(A)Returns vector of column means.colSums(A)Returns vector of coumn sums.2.4.7.其它函数交叉乘积(cross product), 函数为crossprod( )，crossprod(X,Y)表示一般的内积XY，即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中除去(sweep out)。2.5.1 求均值> 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) # 求各列的均值1 0.3386202 0.7320669 -0.4624578 -0.32254602.5.2 标准化> scale(m, center=T, scale=T)2.5.3 减去中位数> row.med <- apply(m, MARGIN=1, FUN=median)> sweep(m, MARGIN=1, STATS=row.med, FUN=”-”)第三章 多元正态分布及参数的估计3.1 绘制二元正态密度函数及其相应等高线图书上例2.2.2，时的二元正态密度函数及其等高线图： x<-seq(-3,3,by=0.1) y<-x f<-function(x,y,a=1,b=1,r=0) a1=sqrt(a) b1=sqrt(b) d=1-r*r d1=sqrt(d)*a1*b1 z=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 多元正态分布的参数估计3.2.1 多元正态总体的相关量设观测数据阵为 样本均值向量设，=1，2，则样本均值向量Xn：,由2.5.1可得：> 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)/n> A<-t(x)%*%m%*%xA即为所求样本离差阵。 样本协方差阵R中求样本协方差阵的函数为cov（）。样本数据阵X的协方差矩阵S即为：> S<-cov(X) 样本相关阵R中求样本协方差阵的函数为cor（）。样本数据阵X的协方差矩阵R即为:> R<-cor(X)3.2.2 极大似然估计极大似然估计法是建立在极大似然原理基础上的一种统计方法。设总体X，其概率密度函数（连续情况）或分布律（离散情况）为，其中是未知参数(或未知参数向量)。设X1，X2，Xn为取自总体X的样本，则似然函数为： ,)=求使似然函数达到最大的参数的值，即极大似然估计值。在单参数场合，在R中可以使用函数optimize( )求极大似然估计值。optimize( )的调用格式如下：optimize(f = , interval = , lower = min(interval), upper = max(interval), maximum = TRUE, tol = .Machine$double.eps0.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 单总体均值向量的检验及置信域4.2.1均值向量的检验书上例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.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-u0> mm<-diag(1,20)-matrix(1,20,20)/n> a<-t(x)%*%mm%*%x # 样本离差阵> ai=solve(a) > dd=xm%*%ai%*%t(xm)> d2=(n-1)*dd> t2=n*d2;> f<-(n-p)*t2/(n-1)*p) # 检验统计量> f ,11, 2.904546> fa<-qf(0.95,p,n-p) # 自由度为(p,n-p)的F分布的0.95分位数> fa1 3.196777> b<-1-pf(f,p,n-p) # 尾概率值> b ,11, 0.06492834> beta<-pf(fa,p,n-p,t2) # 犯第二类错误的概率（假设总体均值）> beta1 0.3616381取检验水平为0.05，由尾概率值p=0.064928340.05=，可得相容；同样由F=2.9045463.196777=Fa，也可得相容。在这种情况下，可能犯第二类错误，概率为=0.3616（假定总体均值）。4.2.2样本协方差阵的特征值和特征向量书上例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 ,31, 2.879368 10.0100 -1.8090532, 10.010000 199.7884 -5.6400003, -1.809053 -5.6400 3.627658> a<-eigen(s)> a$values1 200.462464 4.531591 1.301392$vectors ,1 ,2 ,31, -0.05084144 -0.57370364 0.817483512, -0.99828352 0.05302042 -0.024876553, 0.02907156 0.81734508 0.575414524.3多总体均值向量的检验4.3.1 两正态总体均值向量的检验书上例3.3.1，R程序为：> n<-10> m<-10> p<-4> x<-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%*%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)/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 ,41, 64 43 30.5 63> y0 ,1 ,2 ,3 ,41, 51.5 51 40 70.5> a1 ,1 ,2 ,3 ,41, 490 -170 -120.0 -2452, -170 510 10.0 3103, -120 10 322.5 2604, -245 310 260.0 510> a2 ,1 ,2 ,3 ,41, 502.5 60 175 -7.52, 60.0 390 50 195.03, 175.0 50 450 -100.04, -7.5 195 -100 322.5> d2 ,11, 5.972499> t2 ,11, 29.86250> f ,11, 6.221353> pp ,11, 0.003705807取检验水平为0.01，根据尾概率值p=0.0037058070.01=，可得应否定。4.3.2 多个正态总体均值向量的检验-多元方差分析书上例3.3.2，可利用类似例3.2.1或例3.3.1的程序进行计算得出结论。下面我们用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) Response x1 : Df Sum Sq Mean Sq F value Pr(>F) rate 2 39066 19533 8.878 0.0004401 *Residuals 57 125409 2200 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Response x2 : Df Sum Sq Mean Sq F value Pr(>F) rate 2 4017 2009 2.8293 0.06738 .Residuals 57 40467 710 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Respo
