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