1.4　多元正態分布

正態分布在統計學中有舉足輕重的地位。一方面是因為它可以用來描述很多自然現象，另一方面它又是許多其他理論分布的漸近或者極限。多元正態分布可以看作一元正態分布的推廣。很多多元統計分析的理論都是建立在多元正態總體的基礎上。

1.4.1　多元正態分布的定義

在概率論中，一元正態分布的密度函數為

（1.14）

為了便于推廣，我們將上式改寫為

（1.15）

其中，表示的轉置。這里，因為和都是一維的數，所以與相同。可以看作以標準差為測量單位的到的距離的平方。

將式（1.15）中的一元距離推廣到多元距離，并對前面的系數進行相應的變換，就可以得到多元正態密度函數。

設p元隨機向量X = 的概率密度函數為

（1.16）

則稱X服從p元正態分布，記作，也稱X為p元正態變量。其中，為協方差矩陣的行列式。

MASS包的mvrnorm()函數可以生成服從多元正態分布的隨機數向量。例如，下面的命令生成了一個均值向量為（2, 3, 5）、協方差矩陣為3階單位矩陣、長度為50的服從多元正態分布的隨機數向量。為了讓結果具有重復性，我們設定了生成隨機數的種子。

>mu<-c(2,3,5)#設定均值向量
>sigma<-diag(3)#設定協方差矩陣
>n<-50#設定樣本量
>set.seed(1234)#設定隨機數種子
>library(MASS)
>m<-mvrnorm(n,mu,sigma)
>head(m)
[,1][,2][,3]
[1,]1.2895933.2533194.439524
[2,]2.2568842.9714534.769823
[3,]1.7533082.9571306.558708
[4,]1.6524574.3686025.070508
[5,]1.0483812.7742295.129288
[6,]1.9549724.5164716.715065

1.4.2　多元正態分布的檢驗

要檢驗一個隨機向量是否服從多元正態分布，可以借助MVN包里的函數mvn()。例如，檢驗1.4.1節生成的隨機向量是否服從多元正態分布：

>library(MVN)
>mvn(m)
$multivariateNormality
TestHZpvalueMVN
1Henze-Zirkler0.54212440.7899205YES
　
$univariateNormality
TestVariableStatisticpvalueNormality
1Anderson-DarlingColumn10.31260.5379YES
2Anderson-DarlingColumn20.18980.8948YES
3Anderson-DarlingColumn30.21000.8529YES
　
$Descriptives
nMeanStd.DevMedianMinMax25th
1501.7461000.98933391.695130-0.053247224.1001091.049417
2503.1464080.90544723.1525790.690831125.1873332.638703
3505.0344040.92587004.9273603.033382847.1689564.440683
75thSkewKurtosis
12.2900860.41177454-0.38082462
23.629436-0.043533180.04864663
35.6981770.16269192-0.52369419

函數mvn()中，參數mvnTest的默認值為“hz”，即對隨機向量做Henze-Zirkler多元正態性檢驗；參數univariateTest默認為“AD”，即進行單變量的Anderson-Darling正態分布擬合優度檢驗。根據上面的輸出結果可以認為，各列變量均服從單變量的正態分布，且3個變量服從多元正態分布（P = 0.79）。

需要強調的是，單變量的正態性是多變量多元正態性的必要非充分條件。即使每個變量的分布呈正態分布，所有變量的組合也未必呈多元正態分布。例如，對于表1-2中的4項臨床檢測指標，檢測它們的正態性：

>mvn(bio,univariateTest="SW")
$multivariateNormality
TestHZpvalueMVN
1Henze-Zirkler1.0473870.00628909NO
$univariateNormality
TestVariableStatisticpvalueNormality
1Shapiro-WilkFIB0.95350.1348YES
2Shapiro-WilklnPT0.94800.0903YES
3Shapiro-WilkPTA0.98530.9041YES
4Shapiro-WilklnCHE0.94180.0579YES
$Descriptives
nMeanStd.DevMedianMinMax25th75th
FIB362.4700000.83373862.390.785.102.14752.8325
lnPT362.6780560.19223232.642.303.082.56752.7525
PTA3679.77777823.025589679.0032.00128.0067.750094.7500
lnCHE368.1330560.52223058.197.158.977.67758.5975
SkewKurtosis
FIB0.60465131.2879970
lnPT0.5042816-0.3325948
PTA-0.1565475-0.5762517
lnCHE-0.2425023-1.2648804

在上面的命令中，我們將參數univariateTest設為“SW”，即進行單變量的Shapiro-Wilk正態性檢驗。檢驗結果表明，所有P值均大于0.05，即每個變量都服從正態分布。但是，Henze-Zirkler檢驗的結果表明這4個指標不服從多元正態分布（P = 0.006 < 0.05）。輸出結果還給出了4項指標的常用描述性統計量。

1.4.3　二元正態分布及其參考值范圍

在式（1.16）中，如果p = 2，則稱X服從二元正態分布。設和的均值分別為和，方差分別為和，和的相關系數為，則X的概率密度函數表達式為

（1.17）

為了直觀地展示二元正態分布，我們用下面的命令繪制了均值向量為（0, 0）、方差都為1、相關系數為0的二元正態分布的三維分布曲面。

>mu1<-0;mu2<-0#設定均值
>s1<-1;s2<-1#設定方差
>rho<-0#設定相關系數
>#定義密度函數
>dens<-function(x,y){
+(2*pi*s1*s2*sqrt(1-rho^2))^-1*
+exp(-0.5*(1-rho^2)^-1*
+((x-mu1)^2/s1^2-
+2*rho*(x-mu1)*(y-mu2)/(s1*s2)+
+(y-mu2)^2/s2^2))
+}
>x<-seq(-3,3,length=50)
>y<-seq(-3,3,length=50)
>z<-outer(x,y,dens)#計算z=f(x,y),并且輸出給z
>persp(x,y,z,theta=55,phi=25)#繪制二元正態分布曲面

繪圖結果如圖1-1所示。讀者可以重新設定不同的參數繪制相應的二元正態分布曲面。

圖1-1 二元正態分布曲面（,,）

對于二元正態分布，在某一個變量的取值固定時，另外一個變量服從（單變量的）正態分布。因此，從圖1-1的剖面可以看到很多條正態分布曲線。

從一元正態分布很容易得到變量的95%參考值范圍。一般來說，由于多元正態分布確定的多元參考值范圍考慮了多個指標之間的相關性，因此要比單個指標的參考值范圍的簡單聯合更合理。對于服從二元正態分布的變量，它們的參考值范圍在幾何上是一個橢圓，稱為置信橢圓（confidence ellipse）。car包中的函數dataEllipse()可以用于繪制置信橢圓，例如：

#install.packages("car")
>library(car)
>ellipse<-dataEllipse(cirr$FIB,cirr$PTA,levels=0.95,
+lwd=1.5,center.cex=1,
+xlim=c(0,5),ylim=c(10,140))

繪圖結果如圖1-2（a）所示。該置信橢圓表示，根據樣本數據，肝硬化患者的FIB和PTA的檢測值落在該區域的約占95%。

ggplot2包也可以用于繪制置信橢圓，繪圖結果如圖1-2（b）所示。

#install.packages("ggplot2")
>library(ggplot2)
>ggplot(cirr,aes(FIB,PTA))+geom_point(color='red')+
+stat_ellipse(type='norm')

（a）

（b）

圖1-2 二維相關數據FIB和PTA的95%參考值范圍

使用car包繪制置信橢圓的好處是可以提取橢圓上各點的坐標。查看對象ellipse的前6行，結果如下。

>head(ellipse)
xy
[1,]4.060018138.6389
[2,]4.222379138.1928
[3,]4.358175136.8611
[4,]4.465349134.6641
[5,]4.542276131.6351
[6,]4.587789127.8200