R-多元統(tǒng)計分析上機講義

應用多元統(tǒng)計分析

R實驗上機講義TOC\o"1-5"\h\z應用多元統(tǒng)計分析4AppliedMultivariateStatisticalAnalysis4第一章緒論4第二章矩陣4矩陣的建立4矩陣的下標（index）與子集（元素）的提取6矩陣四則運算7矩陣的加減運算7矩陣的相乘8矩陣的求逆8矩陣的其他一些代數運算8求轉置矩陣8提取對角元素8矩陣的合并與拉直8方陣的行列式9矩陣的特征根和特征向量9其它函數9矩陣的統(tǒng)計運算11求均值11標準化11減去中位數11第三章多元正態(tài)分布及參數的估計12繪制二元正態(tài)密度函數及其相應等高線圖12多元正態(tài)分布的參數估計14多元正態(tài)總體的相關量14極大似然估計14第四章多元正態(tài)總體參數的假設檢驗15幾個重要統(tǒng)計量的分布15單總體均值向量的檢驗及置信域16均值向量的檢驗16樣本協(xié)方差陣的特征值和特征向量17多總體均值向量的檢驗17兩正態(tài)總體均值向量的檢驗17多個正態(tài)總體均值向量的檢驗-多元方差分析19協(xié)方差陣的檢驗204.4.2多總體協(xié)方差陣的檢驗20獨立性檢驗20正態(tài)性檢驗21第五章判別分析22距離判別22馬氏距離22兩總體的距離判別22多個總體的距離判別26貝葉斯判別法及廣義平方距離判別法26先驗用一率（先知知識）26廣義平方距離26后驗概率（條件概率）27貝葉斯判別準則27費希爾（Fisher）判別29第六章聚類分析30距離和相似系數30距離31數據中心化與標準化變換31相似系數31系統(tǒng)聚類法31類個數的確定34動態(tài)聚類法36變量聚類方法36第七章主成分分析37樣本的主成分38主成分分析的應用39第八章因子分析42參數估計方法42方差最大的正交旋轉45因子得分45第九章對應分析方法46第十章典型相關分析48應用多元統(tǒng)計分析AppliedMultivariateStatisticalAnalysis第一章緒論在實際問題中，很多隨機現象涉及到的變量不是一個，而是經常是多個變量，并且這些變量間又存在一定的聯(lián)系。我們經常需要處理多個變量的觀測數據，如果用一元統(tǒng)計方法，由于忽視了各個變量之間可能存在的相關性，一般說來，丟失信息太多，分析的結果不能客觀全面反映數據所包含的內容，因此，我們就需要用到多元統(tǒng)計的方法。多元統(tǒng)計分析(MultivariateStatisticalAnalysis)也稱多變量統(tǒng)計分析、多因素統(tǒng)計分析或多元分析，是研究客觀事物中多變量(多因素或多指標)之間的相互關系和多樣品對象之間差異以及以多個變量為代表的多元隨機變量之間的依賴和差異的現代統(tǒng)計分析理論和方法。多元統(tǒng)計分析是解決實際問題的有效的數據處理方法。隨著電子計算機使用的日益普及，多元統(tǒng)計統(tǒng)計方法已廣泛地應用于自然科學、社會科學的各個方面。第二章矩陣矩陣即是二維的數組，它非常的重要，以至于需要單獨討論。由于矩陣應用非常廣泛，因此對它定義了一些特殊的應用和操作，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"注意：①循環(huán)準則仍然適用于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)與子集(元素)的提取矩陣的下標可以使用正整數、負整數和邏輯表達式，從而實現子集的提取或修改?？疾榫仃噚<-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)不僅低效而且還有一種潛在的不穩(wěn)定性。矩陣的其他一些代數運算求轉置矩陣轉置函數為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()把幾個矩陣列向拼成一個大矩陣，這些矩陣列數應該相同。(如果參與合并的矩陣比其它矩陣行數少或列數少，則循環(huán)不足后合并。)例如：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矩陣的統(tǒng)計運算函數cov()和cor()分別用于計算矩陣的協(xié)方差陣和相關系數陣。矩陣的排列是有方向性的，在R中規(guī)定矩陣是按列排的，若沒有特別說明，函數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=-”)第三章多元正態(tài)分布及參數的估計繪制二元正態(tài)密度函數及其相應等高線圖書上例2.2.2,5=(l￡=Lp=0時的二元正態(tài)密度函數及其等高線圖：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時的二元正態(tài)密度函數")#密度函數圖contour(x,y,z)#等高線圖image(x,y,z)#等高線圖，實際數據大小用不同色彩表示所得圖形為：

-3ct1=o2=1,p=0時的二元正態(tài)密度國數-32-10123-3ct1=o2=1,p=0時的二元正態(tài)密度國數-32-10123相應等高線圖f,把x的任一個元素與f,把x的任一個元素與y的任意一個元素搭配起來作為f的自變量計算得到新的元素值，當函數缺省時表示乘積情況。對參數進行修改，可以繪制任一二元正態(tài)密度函數及其相應的等高線圖。多元正態(tài)分布的參數估計多元正態(tài)總體的相關量設觀測數據陣為亞]…正"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“口為所求樣本離差陣。?樣本協(xié)方差陣R中求樣本協(xié)方差陣的函數為cov()。樣本數據陣X勺協(xié)方差矩陣S即為：>S<-cov(X)?樣本相關陣R中求樣本協(xié)方差陣的函數為cor()。樣本數據陣X勺協(xié)方差矩陣R即為：>R<-cor(X)極大似然估計極大似然估計法是建立在極大似然原理基礎上的一種統(tǒng)計方法。設總體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,…)三者主要區(qū)別是：函數nlm()僅使用牛頓-拉夫遜算法求函數的最小值點；函數0Ptim()提供method選項給出的5種方法中的一種進行優(yōu)化；上面二個可用于多維函數的極值問題，，而函數optimize()僅適用于一維函數，但可以用于最大與最小值點。(具體選項見幫助。)第四章多元正態(tài)總體參數的假設檢驗在一元統(tǒng)計中，用于檢驗一元正態(tài)總體參數N,小3的抽樣分布有好分布，t分布、F分布風，它們都是來自總體NS，口弓的隨機樣本導出的檢驗統(tǒng)計量。推廣到多元正態(tài)總體后，也有相應于以上三個常用分布的統(tǒng)計量：威沙特(Wishart)統(tǒng)計量，霍特林(Hotelling)T?統(tǒng)計量，威爾克斯(Wilks)A統(tǒng)計量，這些統(tǒng)計量是多元統(tǒng)計分析所涉及的假設檢驗問題的基礎。幾個重要統(tǒng)計量的分布對于多元正態(tài)總體來說，存在幾個重要的統(tǒng)計量：威沙特(Wishart)統(tǒng)計量，霍特林（Hotelling）T?統(tǒng)計量，威爾克斯（Wilks）A統(tǒng)計量等，討論這些統(tǒng)計量的分布是多元統(tǒng)計分析所涉及的假設檢驗問題的基礎。單總體均值向量的檢驗及置信域均值向量的檢驗書上例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)#檢驗統(tǒng)計量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)。樣本協(xié)方差陣的特征值和特征向量書上例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多總體均值向量的檢驗兩正態(tài)總體均值向量的檢驗書上例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