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

共享知識(shí) 分享快樂應(yīng)用多元統(tǒng)計(jì)分析R實(shí)驗(yàn)上機(jī)講義卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂應(yīng)用多元統(tǒng)計(jì)分析...........................................................................................................................4AppliedMultivariateStatisticalAnalysis...................................................................................4第一章緒論..............................................................................................................................4第二章矩陣..............................................................................................................................42.1矩陣的建立........................................................................................................................42.2矩陣的下標(biāo)(index)與子集(元素)的提取...........................................................................62.3矩陣四則運(yùn)算...................................................................................................................72.3.1矩陣的加減運(yùn)算.....................................................................................................72.3.2矩陣的相乘............................................................................................................82.3.3矩陣的求逆............................................................................................................82.4矩陣的其他一些代數(shù)運(yùn)算.................................................................................................82.4.1求轉(zhuǎn)置矩陣............................................................................................................82.4.2提取對(duì)角元素.........................................................................................................82.4.3矩陣的合并與拉直..................................................................................................82.4.4方陣的行列式..........................................................................................................92.4.5矩陣的特征根和特征向量.....................................................................................92.4.6其它函數(shù)................................................................................................................92.5矩陣的統(tǒng)計(jì)運(yùn)算..............................................................................................................112.5.1求均值..................................................................................................................112.5.2標(biāo)準(zhǔn)化..................................................................................................................112.5.3減去中位數(shù)..........................................................................................................11第三章多元正態(tài)分布及參數(shù)的估計(jì).....................................................................................123.1繪制二元正態(tài)密度函數(shù)及其相應(yīng)等高線圖..................................................................123.2多元正態(tài)分布的參數(shù)估計(jì)..............................................................................................143.2.1多元正態(tài)總體的相關(guān)量.......................................................................................143.2.2極大似然估計(jì).......................................................................................................14第四章多元正態(tài)總體參數(shù)的假設(shè)檢驗(yàn).................................................................................154.1幾個(gè)重要統(tǒng)計(jì)量的分布..................................................................................................154.2單總體均值向量的檢驗(yàn)及置信域..................................................................................164.2.1均值向量的檢驗(yàn)....................................................................................................164.2.2樣本協(xié)方差陣的特征值和特征向量....................................................................174.3多總體均值向量的檢驗(yàn)...................................................................................................174.3.1兩正態(tài)總體均值向量的檢驗(yàn)...............................................................................174.3.2多個(gè)正態(tài)總體均值向量的檢驗(yàn)-多元方差分析..................................................194.4協(xié)方差陣的檢驗(yàn)...............................................................................................................204.4.2多總體協(xié)方差陣的檢驗(yàn).......................................................................................204.5獨(dú)立性檢驗(yàn)......................................................................................................................204.6正態(tài)性檢驗(yàn)......................................................................................................................21第五章判別分析.....................................................................................................................225.1距離判別..........................................................................................................................225.1.1馬氏距離..............................................................................................................225.1.2兩總體的距離判別...............................................................................................225.1.3多個(gè)總體的距離判別...........................................................................................265.2貝葉斯判別法及廣義平方距離判別法...........................................................................26卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂5.2.1先驗(yàn)概率(先知知識(shí))............................................................................................265.2.2廣義平方距離.......................................................................................................265.2.3后驗(yàn)概率（條件概率）.......................................................................................275.2.4貝葉斯判別準(zhǔn)則...................................................................................................275.3費(fèi)希爾（Fisher）判別.....................................................................................................29第六章聚類分析.....................................................................................................................306.2距離和相似系數(shù)...............................................................................................................306.2.1距離.......................................................................................................................316.2.2數(shù)據(jù)中心化與標(biāo)準(zhǔn)化變換....................................................................................316.2.3相似系數(shù)...............................................................................................................316.3系統(tǒng)聚類法.....................................................................................................................316.4類個(gè)數(shù)的確定..................................................................................................................346.5動(dòng)態(tài)聚類法......................................................................................................................366.7變量聚類方法..................................................................................................................36第七章主成分分析.................................................................................................................377.2樣本的主成分.................................................................................................................387.3主成分分析的應(yīng)用..........................................................................................................39第八章因子分析.....................................................................................................................428.3參數(shù)估計(jì)方法.................................................................................................................428.4方差最大的正交旋轉(zhuǎn)......................................................................................................458.5因子得分.........................................................................................................................45第九章對(duì)應(yīng)分析方法.............................................................................................................46第十章典型相關(guān)分析.............................................................................................................48卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂應(yīng)用多元統(tǒng)計(jì)分析AppliedMultivariateStatisticalAnalysis第一章緒論在實(shí)際問題中，很多隨機(jī)現(xiàn)象涉及到的變量不是一個(gè)， 而是經(jīng)常是多個(gè)變量， 并且這些變量間又存在一定的聯(lián)系。我們經(jīng)常需要處理多個(gè)變量的觀測(cè)數(shù)據(jù)，如果用一元統(tǒng)計(jì)方法，由于忽視了各個(gè)變量之間可能存在的相關(guān)性， 一般說來，丟失信息太多，分析的結(jié)果不能客觀全面反映數(shù)據(jù)所包含的內(nèi)容，因此，我們就需要用到多元統(tǒng)計(jì)的方法。多元統(tǒng)計(jì)分析 (MultivariateStatisticalAnalysis) 也稱多變量統(tǒng)計(jì)分析、多因素統(tǒng)計(jì)分析或多元分析，是研究客觀事物中多變量 (多因素或多指標(biāo) )之間的相互關(guān)系和多樣品對(duì)象之間差異以及以多個(gè)變量為代表的多元隨機(jī)變量之間的依賴和差異的現(xiàn)代統(tǒng)計(jì)分析理論和方法。 多元統(tǒng)計(jì)分析是解決實(shí)際問題的有效的數(shù)據(jù)處理方法。 隨著電子計(jì)算機(jī)使用的日益普及， 多元統(tǒng)計(jì)統(tǒng)計(jì)方法已廣泛地應(yīng)用于自然科學(xué)、社會(huì)科學(xué)的各個(gè)方面。第二章矩陣矩陣即是二維的數(shù)組，它非常的重要，以至于需要單獨(dú)討論。由于矩陣應(yīng)用非常廣泛，因此對(duì)它定義了一些特殊的應(yīng)用和操作，R包括許多只對(duì)矩陣操作的操作符和函數(shù)。2.1矩陣的建立在R中最為常用的是用命令 matrix() 建立矩陣，而對(duì)角矩陣常用函數(shù) 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卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂X<-matrix(1:4,2)#等價(jià)于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"注意：①循環(huán)準(zhǔn)則仍然適用于matrix()，但要求數(shù)據(jù)項(xiàng)的個(gè)數(shù)等于矩陣的列數(shù)的倍數(shù),否則會(huì)出現(xiàn)警告。②矩陣的維數(shù)使用 c()會(huì)得到不同的結(jié)果 (除非是方陣),因此需要小心。③數(shù)據(jù)項(xiàng)填充矩陣的方向可通過參數(shù) byrow來指定, 其缺省是按列填充的byrow=FALSE)，byrow=TRUE表示按行填充數(shù)據(jù)。再看幾個(gè)例子：>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):數(shù)據(jù)長度[4]不是矩陣列數(shù)[3]的整倍數(shù)>X<-matrix(1:4,c(2,3))#不經(jīng)常使用>X[,1][,2][1,]13[2,]24>X<-matrix(1:4,2,4,byrow=TRUE)#按行填充>X[,1][,2][,3][,4]卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí)分享快樂[1,]1234[2,]1234因?yàn)榫仃囀菙?shù)組的特例,R中數(shù)組由函數(shù)array()建立,因此矩陣也可以用函數(shù)array()來建立，其一般格式為：>array(data,dim,dimnames)其中data為一向量，其元素用于構(gòu)建數(shù)組；dim為數(shù)組的維數(shù)向量（為數(shù)值型向量）；dimnames為由各維的名稱構(gòu)成的向量（為字符型向量)，缺省為空。看幾個(gè)例子：>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矩陣的下標(biāo)(index)與子集(元素)的提取矩陣的下標(biāo)可以使用正整數(shù)、 負(fù)整數(shù)和邏輯表達(dá)式， 從而實(shí)現(xiàn)子集的提取或修改。 考查矩陣x<-matrix(1:6,2,3)x[,1][,2][,3][1,] 1 3 5[2,] 2 4 6提取一個(gè)元素x[2,2][1]4提取若一個(gè)或若干個(gè)行或列x[2,2]4>x[2,]246卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂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去掉某一個(gè)或若干個(gè)行與列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矩陣四則運(yùn)算矩陣也可以進(jìn)行四則運(yùn)算(“+”、“-”、“*”、“/”，“^”)，分別解釋為矩陣對(duì)應(yīng)元素的四則運(yùn)算。在實(shí)際應(yīng)用中，比較有實(shí)際應(yīng)用的是矩陣的相加，相減，相乘和矩陣的求逆。矩陣的加減運(yùn)算一般要求矩陣形狀完全相同（dim屬性完全相同），矩陣的相乘一般要求一矩陣的列維數(shù)與另一矩陣的行維數(shù)相同，而矩陣要求逆的話，一般要求它為一方陣。 矩陣的加減運(yùn)算若A，B為兩個(gè)形狀相同的矩陣，兩矩陣的和為 C，R中表達(dá)式為：C<-A+B兩矩陣的差為 D，R中表達(dá)式為：D<-A-B卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂矩陣也可以與數(shù)進(jìn)行加減， A+5表示A中的每個(gè)元素加上 5。 矩陣的相乘操作符%*% 用于矩陣相乘。若矩陣 A的列數(shù)等于矩陣 B的行數(shù)，矩陣 A乘以矩陣表示為：A%*%B注：X*Y表示兩個(gè)矩陣的逐元相乘，而不是 X和Y的乘積。 矩陣的求逆若矩陣A為一方陣，矩陣的逆可以用下面的命令計(jì)算： solve(A)。操作符solve()可以用來求解線性方程組： Ax=b，解為solve(A,b)在數(shù)學(xué)上，用直接求逆的辦法解x<-solve(A)%*%b相比solve(A,b)不僅低效而且還有一種潛在的不穩(wěn)定性。2.4矩陣的其他一些代數(shù)運(yùn)算 求轉(zhuǎn)置矩陣轉(zhuǎn)置函數(shù)為 t() ，矩陣X的轉(zhuǎn)置為 t(X)。 提取對(duì)角元素提取對(duì)角元的函數(shù)為 diag()。例如：X<-matrix(1:4,2,2)diag(X)[1]14事實(shí)上，diag()的作用依賴于自變量， diag（vector）返回以自變量（向量）為主對(duì)角元素的對(duì)角矩陣；diag（matrix）返回由矩陣的主對(duì)角元素所組成的向量；diag（k）（k為標(biāo)量）返回k階單位陣。矩陣的合并與拉直函數(shù)cbind（）把幾個(gè)矩陣橫向拼成一個(gè)大矩陣，這些矩陣行數(shù)應(yīng)該相同；函數(shù) rbind（）把幾個(gè)矩陣列向拼成一個(gè)大矩陣，這些矩陣列數(shù)應(yīng)該相同。 （如果參與合并的矩陣比其它矩陣行數(shù)少或列數(shù)少，則循環(huán)不足后合并。）例如：>m1<-matrix(1,nr=2,nc=2)>m1[,1][,2][1,] 1 1[2,] 1 1m2<-matrix(2,nr=2,nc=2)m2卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂[,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 矩陣的特征根和特征向量函數(shù)eigen()用來計(jì)算矩陣的特征值和特征向量。這個(gè)函數(shù)的返回值是一個(gè)含有values和vectors兩個(gè)分量的列表。命令A(yù)<-eigen(X)>A$values[1]5.3722813-0.3722813$vectors[,1][,2][1,]-0.5657675-0.9093767[2,]-0.82456480.4159736Inthefollowingexamples,AandBarematricesandxandbareavectors.OperatororFunctionDescriptionA*B卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂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,...)卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂Combinematrices(vectors)vertically.Returnsamatrix.rowMeans(A)Returnsvectorofrowmeans.rowSums(A)Returnsvectorofrowsums.colMeans(A)Returnsvectorofcolumnmeans.colSums(A)Returnsvectorofcoumnsums.其它函數(shù)交叉乘積(crossproduct), 函數(shù)為 crossprod() ，crossprod(X,Y) 表示一般的內(nèi)積X′Y，即X的每一列與Y的每一列的內(nèi)積組成的矩陣；QR分解,函數(shù)為qr(),矩陣X的QR分解為X=QR，Q為正交陣，R為上三角陣；等等。2.5矩陣的統(tǒng)計(jì)運(yùn)算函數(shù)cov()和cor()分別用于計(jì)算矩陣的協(xié)方差陣和相關(guān)系數(shù)陣。矩陣的排列是有方向性的，在R中規(guī)定矩陣是按列排的，若沒有特別說明，函數(shù)max()，min()，median()，var()，sd()，sum()，cumsum()，cumprod()，cummax()，cummin()的使用對(duì)于矩陣也是按列計(jì)算的，但也可以通過選項(xiàng)MARGIN來改變。下面我們要用到對(duì)一個(gè)對(duì)象施加某種運(yùn)算的函數(shù)apply()，其格式為>apply(X,MARGIN,FUN)其中X為參與運(yùn)算的矩陣 ,FUN為上面的一個(gè)函數(shù)或“ +”、“-”、“*”、“＼”(必須放在引號(hào)中)，MARGIN=1表示按列計(jì)算， MARGIN=2表示按行計(jì)算。我們還用到sweep()函數(shù)，命令>sweep(X,MARGIN,STATS,FUN)表示從矩陣X中按MATGIN計(jì)算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 標(biāo)準(zhǔn)化>scale(m,center=T,scale=T) 減去中位數(shù)>row.med<-apply(m,MARGIN=1,FUN=median)>sweep(m,MARGIN=1,STATS=row.med,FUN= -”)卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂第三章多元正態(tài)分布及參數(shù)的估計(jì)3.1繪制二元正態(tài)密度函數(shù)及其相應(yīng)等高線圖書上例， 時(shí)的二元正態(tài)密度函數(shù)及其等高線圖：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) #外積函數(shù)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時(shí)的二元正態(tài)密度函數(shù) ")密度函數(shù)圖contour(x,y,z) # 等高線圖image(x,y,z) # 等高線圖，實(shí)際數(shù)據(jù)大小用不同色彩表示所得圖形為：卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂相應(yīng)等高線圖Outer(x,y,f)是一個(gè)一般性的外積函數(shù)，調(diào)用函數(shù)f，把x的任一個(gè)元素與y的任意一個(gè)元素搭配起來作為f的自變量計(jì)算得到新的元素值，當(dāng)函數(shù)缺省時(shí)表示乘積情況。卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂對(duì)參數(shù)進(jìn)行修改，可以繪制任一二元正態(tài)密度函數(shù)及其相應(yīng)的等高線圖。3.2多元正態(tài)分布的參數(shù)估計(jì) 多元正態(tài)總體的相關(guān)量設(shè)觀測(cè)數(shù)據(jù)陣為樣本均值向量設(shè) ，=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即為所求樣本離差陣。樣本協(xié)方差陣R中求樣本協(xié)方差陣的函數(shù)為 cov（）。樣本數(shù)據(jù)陣 X的協(xié)方差矩陣S即為：>S<-cov(X)樣本相關(guān)陣R中求樣本協(xié)方差陣的函數(shù)為 cor（）。樣本數(shù)據(jù)陣 X的協(xié)方差矩陣R即為:>R<-cor(X) 極大似然估計(jì)極大似然估計(jì)法是建立在極大似然原理基礎(chǔ)上的一種統(tǒng)計(jì)方法。設(shè)總體 X，其概率密度函數(shù)（連續(xù)情況）或分布律（離散情況）為 ，其中 是未知參數(shù)(或未知參數(shù)向量 )。設(shè)X1，X2，?，Xn為取自總體X的樣本，則似然函數(shù) 為：卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂?, )=求使似然函數(shù)達(dá)到最大的參數(shù) 的值，即極大似然估計(jì)值。在單參數(shù)場合，在R中可以使用函數(shù)optimize()求極大似然估計(jì)值。optimize()的調(diào)用格式如下：optimize(f=,interval=,lower=min(interval),upper=max(interval),maximum=TRUE,tol= .Machine$double.eps^0.25, ?)說明：f是似然函數(shù)， interval是參數(shù) 的取值范圍， lower是 的下界，upper是 的上界，maximum=TRUE是求極大值，否則(maximum=FALSE)表示求函數(shù)的極小值，tol是表示求值的精確度，?是對(duì)f的附加說明。在多參數(shù)場合，在R中用函數(shù)optim()或者nlm()來求似然函數(shù)的極大值，并求相應(yīng)的極大值點(diǎn)。optim()的調(diào)用格式如下：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, ?)三者主要區(qū)別是：函數(shù) nlm()僅使用牛頓-拉夫遜算法求函數(shù)的最小值點(diǎn)；函數(shù) optim()提供method選項(xiàng)給出的5種方法中的一種進(jìn)行優(yōu)化；上面二個(gè)可用于多維函數(shù)的極值問題 ,，而函數(shù)optimize()僅適用于一維函數(shù)，但可以用于最大與最小值點(diǎn)。（具體選項(xiàng)見幫助。）第四章多元正態(tài)總體參數(shù)的假設(shè)檢驗(yàn)在一元統(tǒng)計(jì)中，用于檢驗(yàn)一元正態(tài)總體參數(shù) ， 的抽樣分布有 分布， 分布、F分布風(fēng)，它們都是來自總體 的隨機(jī)樣本導(dǎo)出的檢驗(yàn)統(tǒng)計(jì)量。推廣到多元正態(tài)總體后，也有相應(yīng)于以上三個(gè)常用分布的統(tǒng)計(jì)量：威沙特（ Wishart）統(tǒng)計(jì)量，霍特林（ Hotelling ）統(tǒng)計(jì)量，威爾克斯（ Wilks） 統(tǒng)計(jì)量，這些統(tǒng)計(jì)量是多元統(tǒng)計(jì)分析所涉及的假設(shè)檢驗(yàn)問題的基礎(chǔ)。4.1幾個(gè)重要統(tǒng)計(jì)量的分布對(duì)于多元正態(tài)總體來說 ,存在幾個(gè)重要的統(tǒng)計(jì)量 : 威沙特（Wishart）統(tǒng)計(jì)量，霍特林卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂（Hotelling ） 統(tǒng)計(jì)量，威爾克斯（ Wilks） 統(tǒng)計(jì)量等,討論這些統(tǒng)計(jì)量的分布是多元統(tǒng)計(jì)分析所涉及的假設(shè)檢驗(yàn)問題的基礎(chǔ)。4.2單總體均值向量的檢驗(yàn)及置信域均值向量的檢驗(yàn)書上例，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)# 檢驗(yàn)統(tǒng)計(jì)量>f[,1][1,]2.904546>fa<-qf(0.95,p,n-p)# 自由度為(p,n-p) 的F分布的0.95分位數(shù)>fa[1]3.196777>b<-1-pf(f,p,n-p)# 尾概率值>b[,1][1,]0.06492834>beta<-pf(fa,p,n-p,t2)# 犯第二類錯(cuò)誤的概率（假設(shè)總體均值 ）>beta[1]0.3616381取檢驗(yàn)水平為 0.05，由尾概率值 p=0.06492834 0.05= ，可得 相容；同樣由卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂F=2.904546 3.196777=Fa，也可得 相容。在這種情況下，可能犯第二類錯(cuò)誤，概率為=0.3616（假定總體均值 ）。樣本協(xié)方差陣的特征值和特征向量書上例，R程序?yàn)椋簒<-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多總體均值向量的檢驗(yàn) 兩正態(tài)總體均值向量的檢驗(yàn)書上例，R程序?yà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卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂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]卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂[1,]0.003705807取檢驗(yàn)水平為 0.01，根據(jù)尾概率值 p=0.003705807 0.01= ，可得應(yīng)否定 。 多個(gè)正態(tài)總體均值向量的檢驗(yàn) -多元方差分析書上例，可利用類似例 或例 的程序進(jìn)行計(jì)算得出結(jié)論。下面我們用R自帶的manova()函數(shù)進(jìn)行分析。程序如下：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)# 對(duì)每一個(gè)變量進(jìn)行單因素方差分析summary(fit,test="Wilks")# 使用威爾克斯 統(tǒng)計(jì)量程序結(jié)果：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卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂>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結(jié)果說明：（1） 取檢驗(yàn)水平為 0.01，則對(duì)四個(gè)指標(biāo)逐項(xiàng)用一元方差分析方法進(jìn)行檢驗(yàn)，由 p值可得三個(gè)組指標(biāo)間只有第一個(gè)指標(biāo) 有顯著差異( =0.0004401)；（2） 取檢驗(yàn)水平為 0.01，利用威爾克斯 統(tǒng)計(jì)量得到 p=0.003538 0.01，故拒絕原假設(shè)，即認(rèn)為三個(gè)組的指標(biāo)之間有顯著差異。4.4協(xié)方差陣的檢驗(yàn) 多總體協(xié)方差陣的檢驗(yàn)書上例3.4.1，R程序略（類似例3.2.1或例3.3.1）4.5獨(dú)立性檢驗(yàn)書中例，R程序?yàn)椋簒<-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卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂[1]0.02076288取檢驗(yàn)水平為 0.05，根據(jù)尾概率值 p=0.02076288 0.05= ，可得應(yīng)否定 ，由R軟件所的結(jié)果與 SAS軟件所的結(jié)果一致。4.6正態(tài)性檢驗(yàn)書中例，R程序?yàn)椋簒<-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卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂第五章判別分析判別分析是用于判斷樣品所屬類型的一種統(tǒng)計(jì)分析方法。判別分析的目的是對(duì)已知?dú)w類的數(shù)據(jù)建立由數(shù)值指標(biāo)構(gòu)成的歸類規(guī)則，然后把這樣的規(guī)則應(yīng)用到未知?dú)w類的樣品去歸類。在生產(chǎn)、科研和日常生活中經(jīng)常會(huì)遇到如何根據(jù)觀測(cè)到的數(shù)據(jù)資料對(duì)所研究的對(duì)象進(jìn)行判別歸類的問題。判別分析問題一般可以如下描述：設(shè)有k個(gè)維總體，其分布特征已知（如已知分布函數(shù)分別為，，或知道來自各個(gè)總體的訓(xùn)練樣本）。對(duì)給定的一個(gè)新樣品X，判斷它來自哪個(gè)總體。通常我們先對(duì)預(yù)先得到的來自這k個(gè)總體的若干個(gè)樣品(稱為訓(xùn)練樣品)進(jìn)行檢驗(yàn)和歸類，來決定相應(yīng)的判別歸類問題是否有意義及誤判可能性大小。然后再對(duì)給定的一個(gè)或幾個(gè)新的樣品，進(jìn)行判別歸類，即決定它(們)自哪個(gè)總體。解決這個(gè)問題可以有多種途徑,下面我們分別討論幾種常用的方法，如距離判別、貝葉斯判別、Fisher判別等。R通用程序：首先我們要用命令>library(MASS)MASS包里的lda()針對(duì)線性判別分析。加載MASS宏包，再用函數(shù)lda()就可完成判別分析，其基本調(diào)用格式如下：lda(formula,data,...,subset,na.action)說明：formula用法為groupsx1+x2+?，group表明總體來源，x1，x2，?表示分類~指標(biāo)；subset指明訓(xùn)練樣本。具體說明見R幫助。5.1距離判別 馬氏距離馬氏距離定義：樣本X和總體其中 為總體均值向量， 為總體協(xié)方差陣。 兩總體的距離判別判別準(zhǔn)則：其中 = ， 為X到總體的距離。關(guān)于兩總體距離判別的 R程序(參考薛毅教授的《統(tǒng)計(jì)建模與 R軟件》一書)：卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂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表示訓(xùn)練樣本X1，X2，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入），輸入變量TstX是待測(cè)樣本，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入），或向量（一個(gè)待測(cè)樣本）。如果不輸入 TstX（缺省值），則待測(cè)樣本為兩個(gè)訓(xùn)練樣本之和，即訓(xùn)練樣本的回代情況。輸入變量var.equal是邏輯變量，var.equal==TRUE表示兩個(gè)總體的協(xié)方差相同；否則（缺省值）為不同。在上述程序中，用到馬氏距離函數(shù)mahalanobis()，該函數(shù)的使用格式為mahalanobis(x,center,cov,inverted=FALSE,...)其中x是樣本數(shù)據(jù)構(gòu)成的向量或矩陣（ p維），center為樣本中心， cov為樣本的協(xié)方差陣。對(duì)于書中例，調(diào)用discriminiant.distance()進(jìn)行判別（假設(shè)協(xié)方差陣相等）：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),卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂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由程序結(jié)果可得待判樣品 2，3，6，7，8屬于含鉀鹽泉(A盆地)，其余三個(gè)屬于不含鉀鹽泉(B盆地)。利用R自帶函數(shù)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內(nèi)置函數(shù)，可以將 lda() 的輸出應(yīng)用于訓(xùn)練樣品數(shù)據(jù)進(jìn)行預(yù)測(cè)，從而進(jìn)行對(duì)比。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程序結(jié)果:Call:卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂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卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂2.2416596-3.7639282-9.8136273-8.0017623結(jié)果說明：Groupmeans:包含了每組的平均向量；Coefficientsoflineardiscriminants:線性判別系數(shù)；列聯(lián)表表明將訓(xùn)練樣品數(shù)據(jù)代入線性判別函數(shù)后的判別結(jié)果，兩組都沒有錯(cuò)判；由$class可以看出8個(gè)待判樣品，待判樣品2，3，6，7，8屬于含鉀鹽泉(A盆地)，其余三個(gè)屬于不含鉀鹽泉(B盆地)（與上一程序結(jié)果一致）；$posterior給出了后驗(yàn)概率值（具體概念見5.2節(jié)）；6)$x給出了線性判別函數(shù)的數(shù)值。 多個(gè)總體的距離判別類似與兩個(gè)總體的情況，多個(gè)總體的情況，按照距離最近的原則對(duì) X進(jìn)行判別歸類時(shí)，首先計(jì)算樣品到各類的馬氏 (Mahalanobis)距離，然后進(jìn)行比較，把待判樣品判歸距離最小的那個(gè)總體。（自編關(guān)于多個(gè)總體距離判別的 R函數(shù)可參考《統(tǒng)計(jì)建模與 R軟件一書》）。5.2貝葉斯判別法及廣義平方距離判別法 先驗(yàn)概率(先知知識(shí))設(shè)有k個(gè)總體 ，假設(shè)事先對(duì)所研究的問題有一定的認(rèn)識(shí)，這種認(rèn)識(shí)常用先驗(yàn)概率來描述， 即已知這 k個(gè)總體各自出現(xiàn)的概率 （驗(yàn)前概率）為 （顯然 ， =1），這組驗(yàn)前概率 稱為先驗(yàn)概率。 廣義平方距離樣品X到總體 (=1, ?,k)的廣義平方距離 為：,其中是樣品X到總體 的馬氏距離；其中 為第類的組內(nèi)樣本協(xié)方差陣。卑微如螻蟻、堅(jiān)強(qiáng)似大象共享知識(shí) 分享快樂 后驗(yàn)概率（條件概率）當(dāng)樣品X已知時(shí)，它屬于 的概率就稱為后驗(yàn)概率，一般記為 (或 )。 貝葉斯判別準(zhǔn)則幾個(gè)概念：1.錯(cuò)判概率和錯(cuò)判損失；2.關(guān)于先驗(yàn)概率的平均損失。定義5.2.1：設(shè)有k個(gè)總體：，相應(yīng)的先驗(yàn)概率為（，=1）。如果有判別法 ，使得 帶來的平均損失 達(dá)最小，則稱判別法 符合貝葉斯判別準(zhǔn)則，或稱 為貝葉斯判別的解。出于例題需要，學(xué)習(xí)多總體的 Bayes判別程序（兩總體情況參考《統(tǒng)計(jì)建模與 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<-