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

多元分析R實驗上機講義應(yīng)用多元統(tǒng)計分析R實驗上機講義應(yīng)用多元統(tǒng)計分析4Applied Multivariate Statistical Analysis4第一章緒論4第二章矩陣42.1矩陣的建立42.2矩陣的下標(index)與子集(元素)的提取62.3 矩陣四則運算72.3.1 矩陣的加減運算72.3.2 矩陣的相乘82.3.3 矩陣的求逆82.4矩陣的其他一些代數(shù)運算82.4.1 求轉(zhuǎn)置矩陣82.4.2 提取對角元素82.4.3矩陣的合并與拉直82.4.4方陣的行列式92.4.5 矩陣的特征根和特征向量92.4.6 其它函數(shù)92.5 矩陣的統(tǒng)計運算102.5.1 求均值102.5.2 標準化102.5.3 減去中位數(shù)10第三章多元正態(tài)分布及參數(shù)的估計103.1 繪制二元正態(tài)密度函數(shù)及其相應(yīng)等高線圖103.2 多元正態(tài)分布的參數(shù)估計123.2.1 多元正態(tài)總體的相關(guān)量123.2.2 極大似然估計13第四章多元正態(tài)總體參數(shù)的假設(shè)檢驗144.1 幾個重要統(tǒng)計量的分布144.2 單總體均值向量的檢驗及置信域144.2.1均值向量的檢驗144.2.2樣本協(xié)方差陣的特征值和特征向量154.3多總體均值向量的檢驗164.3.1 兩正態(tài)總體均值向量的檢驗164.3.2 多個正態(tài)總體均值向量的檢驗-多元方差分析174.4協(xié)方差陣的檢驗184.4.2 多總體協(xié)方差陣的檢驗184.5獨立性檢驗184.6正態(tài)性檢驗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距離和相似系數(shù)296.2.1距離296.2.2數(shù)據(jù)中心化與標準化變換296.2.3相似系數(shù)306.3 系統(tǒng)聚類法306.4類個數(shù)的確定326.5動態(tài)聚類法346.7變量聚類方法34第七章主成分分析357.2 樣本的主成分367.3 主成分分析的應(yīng)用37第八章因子分析408.3 參數(shù)估計方法408.4 方差最大的正交旋轉(zhuǎn)438.5 因子得分43第九章對應(yīng)分析方法44第十章典型相關(guān)分析46應(yīng)用多元統(tǒng)計分析Applied Multivariate Statistical Analysis第一章 緒論在實際問題中，很多隨機現(xiàn)象涉及到的變量不是一個，而是經(jīng)常是多個變量，并且這些變量間又存在一定的聯(lián)系。我們經(jīng)常需要處理多個變量的觀測數(shù)據(jù)，如果用一元統(tǒng)計方法，由于忽視了各個變量之間可能存在的相關(guān)性，一般說來，丟失信息太多，分析的結(jié)果不能客觀全面反映數(shù)據(jù)所包含的內(nèi)容，因此，我們就需要用到多元統(tǒng)計的方法。多元統(tǒng)計分析(Multivariate Statistical Analysis)也稱多變量統(tǒng)計分析、多因素統(tǒng)計分析或多元分析，是研究客觀事物中多變量(多因素或多指標)之間的相互關(guān)系和多樣品對象之間差異以及以多個變量為代表的多元隨機變量之間的依賴和差異的現(xiàn)代統(tǒng)計分析理論和方法。多元統(tǒng)計分析是解決實際問題的有效的數(shù)據(jù)處理方法。隨著電子計算機使用的日益普及，多元統(tǒng)計統(tǒng)計方法已廣泛地應(yīng)用于自然科學(xué)、社會科學(xué)的各個方面。第二章 矩陣 矩陣即是二維的數(shù)組，它非常的重要，以至于需要單獨討論。由于矩陣應(yīng)用非常廣泛，因此對它定義了一些特殊的應(yīng)用和操作，R 包括許多只對矩陣操作的操作符和函數(shù)。2.1矩陣的建立在R中最為常用的是用命令matrix( )建立矩陣，而對角矩陣常用函數(shù)diag( )建立。例如 X X ,1 ,21, 1 12, 1 1 X 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 X ,1 ,21, 1 32, 2 4 rownames(X) colnames(X) X c da 1 3b 2 4 dim(X)1 2 2 dimnames(X)11 a b21 c d注意：循環(huán)準則仍然適用于matrix( )，但要求數(shù)據(jù)項的個數(shù)等于矩陣的列數(shù)的倍數(shù), 否則會出現(xiàn)警告。矩陣的維數(shù)使用c( )會得到不同的結(jié)果(除非是方陣), 因此需要小心。數(shù)據(jù)項填充矩陣的方向可通過參數(shù)byrow來指定, 其缺省是按列填充的（byrow=FALSE)， byrow=TRUE表示按行填充數(shù)據(jù)。再看幾個例子： X X ,1 ,2 ,3 ,41, 1 3 1 32, 2 4 2 4 X X X ,1 ,21, 1 32, 2 4 X X ,1 ,2 ,3 ,41, 1 2 3 42, 1 2 3 4 因為矩陣是數(shù)組的特例,R中數(shù)組由函數(shù)array( )建立, 因此矩陣也可以用函數(shù)array( )來建立，其一般格式為： array(data, dim, dimnames)其中data為一向量，其元素用于構(gòu)建數(shù)組；dim為數(shù)組的維數(shù)向量（為數(shù)值型向量）；dimnames為由各維的名稱構(gòu)成的向量（為字符型向量)， 缺省為空?？磶讉€例子： A A ,1 ,2 ,31, 1 3 52, 2 4 6 A A ,1 ,2 ,31, 1 3 12, 2 4 2 A A ,1 ,2 ,31, 1 3 52, 2 4 62.2矩陣的下標(index)與子集(元素)的提取矩陣的下標可以使用正整數(shù)、負整數(shù)和邏輯表達式，從而實現(xiàn)子集的提取或修改?？疾榫仃?x 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 x ,1 ,2 ,31, 1 3 NA2, 2 4 NA xis.na(x) x ,1 ,2 ,31, 1 3 12, 2 4 12.3 矩陣四則運算矩陣也可以進行四則運算(“+”、“-”、“*”、“/”，“”)，分別解釋為矩陣對應(yīng)元素的四則運算。在實際應(yīng)用中，比較有實際應(yīng)用的是矩陣的相加，相減，相乘和矩陣的求逆。矩陣的加減運算一般要求矩陣形狀完全相同（dim屬性完全相同），矩陣的相乘一般要求一矩陣的列維數(shù)與另一矩陣的行維數(shù)相同，而矩陣要求逆的話，一般要求它為一方陣。2.3.1 矩陣的加減運算若A，B為兩個形狀相同的矩陣，兩矩陣的和為C，R中表達式為：C-A+B兩矩陣的差為D，R中表達式為：D-A-B矩陣也可以與數(shù)進行加減，A+5表示A中的每個元素加上5。2.3.2 矩陣的相乘操作符%*% 用于矩陣相乘。若矩陣A的列數(shù)等于矩陣B的行數(shù)，矩陣A乘以矩陣B表示為：A%*%B注：X*Y表示兩個矩陣的逐元相乘，而不是X和Y的乘積。2.3.3 矩陣的求逆若矩陣A為一方陣，矩陣的逆可以用下面的命令計算：solve(A)。操作符solve( )可以用來求解線性方程組：Ax=b，解為solve(A,b)在數(shù)學(xué)上，用直接求逆的辦法解x X diag(X)1 1 4事實上，diag( )的作用依賴于自變量，diag（vector）返回以自變量（向量）為主對角元素的對角矩陣；diag（matrix）返回由矩陣的主對角元素所組成的向量；diag（k）（k為標量）返回k階單位陣。2.4.3矩陣的合并與拉直函數(shù)cbind（）把幾個矩陣橫向拼成一個大矩陣，這些矩陣行數(shù)應(yīng)該相同；函數(shù)rbind（）把幾個矩陣列向拼成一個大矩陣，這些矩陣列數(shù)應(yīng)該相同。（如果參與合并的矩陣比其它矩陣行數(shù)少或列數(shù)少，則循環(huán)不足后合并。）例如： m1 m1 ,1 ,21, 1 12, 1 1 m2 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 X ,1 ,21, 1 32, 2 4 det(X)1 -22.4.5 矩陣的特征根和特征向量 函數(shù)eigen( ) 用來計算矩陣的特征值和特征向量。這個函數(shù)的返回值是一個含有values 和vectors 兩個分量的列表。命令A(yù) 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. ABcrossprod(A,B)crossprod(A)AB and AA 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 RR = A.y apply(X, MARGIN, FUN)其中X為參與運算的矩陣, FUN為上面的一個函數(shù)或“+”、“-”、“*”、“”(必須放在引號中)，MARGIN=1表示按列計算，MARGIN=2表示按行計算。我們還用到sweep( )函數(shù)，命令 sweep(X, MARGIN, STATS, FUN)表示從矩陣X中按MATGIN計算STATS，并從X中除去(sweep out)。2.5.1 求均值 m 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 減去中位數(shù) row.med sweep(m, MARGIN=1, STATS=row.med, FUN=”-”)第三章 多元正態(tài)分布及參數(shù)的估計3.1 繪制二元正態(tài)密度函數(shù)及其相應(yīng)等高線圖書上例2.2.2，1=2=1，=0時的二元正態(tài)密度函數(shù)及其等高線圖： 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 Xn ln Xn A m A S R x n p u0 ln x0 xm mm a ai=solve(a) dd=xm%*%ai%*%t(xm) d2=(n-1)*dd t2=n*d2; f f ,11, 2.904546 fa fa1 3.196777 b b ,11, 0.06492834 beta beta1 0.3616381取檢驗水平為=0.05，由尾概率值p=0.064928340.05=，可得H0相容；同樣由F=2.904546 x s s ,1 ,2 ,31, 2.879368 10.0100 -1.8090532, 10.010000 199.7884 -5.6400003, -1.809053 -5.6400 3.627658 a 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 兩正態(tài)總體均值向量的檢驗書上例3.3.1，R程序為： n m p x ln x0 mx a1 y y0 my a2 a xy ai dd d2 t2 f pp 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，根據(jù)尾概率值p=0.0037058070.01=，可得應(yīng)否定H0。4.3.2 多個正態(tài)總體均值向量的檢驗-多元方差分析書上例3.3.2，可利用類似例3.2.1或例3.3.1的程序進行計算得出結(jié)論。下面我們用R自帶的manova()函數(shù)進行分析。程序如下：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 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 Response x3 : Df Sum Sq Mean Sq F value Pr(F)rate 2 13.43 6.72 0.1838 0.8326Residuals 57 2082.50 36.54 Response x4 : Df Sum Sq Mean Sq F value Pr(F)rate 2 17.20 8.60 0.4785 0.6222Residuals 57 1024.40 17.97 summary(fit, test=Wilks) Df Wilks approx F num Df den Df Pr(F) rate 2 0.66212 3.09069 8 108 0.003538 *Residuals 57 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1結(jié)果說明：（1） 取檢驗水平為0.01，則對四個指標逐項用一元方差分析方法進行檢驗，由p值可得三個組指標間只有第一個指標X1有顯著差異(p1=0.0004401)；（2） 取檢驗水平為0.01，利用威爾克斯 統(tǒng)計量得到p=0.003538 x n p x0 xm mm a a0 a1 v b df kc p0 kc1 9.755514 p01 0.02076288取檢驗水平為=0.05，根據(jù)尾概率值p=0.02076288 x n p ln x0 s si m for(i in 1:n)+ xx0-xi,-x0+ dd-xx0%*%si%*%t(xx0)+ print(c(i,dd)+ if(dd=1.386) m m1 7 pp pp1 0.5833333第五章 判別分析判別分析是用于判斷樣品所屬類型的一種統(tǒng)計分析方法。判別分析的目的是對已知歸類的數(shù)據(jù)建立由數(shù)值指標構(gòu)成的歸類規(guī)則，然后把這樣的規(guī)則應(yīng)用到未知歸類的樣品去歸類。 在生產(chǎn)、科研和日常生活中經(jīng)常會遇到如何根據(jù)觀測到的數(shù)據(jù)資料對所研究的對象進行判別歸類的問題。判別分析問題一般可以如下描述：設(shè)有k個m維總體G1，G2，Gk，其分布特征已知（如已知分布函數(shù)分別為F1(x)，F(xiàn)2(x) ，F(xiàn)k(x)，或知道來自各個總體的訓(xùn)練樣本）。對給定的一個新樣品X，判斷它來自哪個總體。通常我們先對預(yù)先得到的來自這k個總體的若干個樣品(稱為訓(xùn)練樣品)進行檢驗和歸類， 來決定相應(yīng)的判別歸類問題是否有意義及誤判可能性大小。然后再對給定的一個或幾個新的樣品，進行判別歸類，即決定它(們)自哪個總體。解決這個問題可以有多種途徑, 下面我們分別討論幾種常用的方法，如距離判別、貝葉斯判別、Fisher判別等。R通用程序：首先我們要用命令library(MASS)MASS包里的lda( )針對線性判別分析。加載MASS宏包，再用函數(shù)lda( )就可完成判別分析，其基本調(diào)用格式如下：lda(formula, data, . , subset, na.action)說明： formula用法為groupsx1 + x2 +，group表明總體來源，x1，x2，表示分類指標；subset指明訓(xùn)練樣本。具體說明見R幫助。5.1距離判別5.1.1 馬氏距離馬氏距離定義：樣本X和總體G 的馬氏距離為：d2X,G=(X-)-1(X-)其中 為總體均值向量， 為總體協(xié)方差陣。5.1.2 兩總體的距離判別判別準則：判 XG1，當(dāng)WX0時，判 XG2，當(dāng)WX0時，待判， 當(dāng)WX0時.其中WX=d2X,G2-d2X,G1，d2X,G為X到總體的距離。關(guān)于兩總體距離判別的R程序(參考薛毅教授的統(tǒng)計建模與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)else if (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