多元正态分布

转自个人微信公众号【Memo_Cleon】的统计学习笔记：多元正态分布检验的R实现方法。

多元正态分布也称多元高斯分布。如同正态分布在单变量分析中的地位类似，在对多个因变量（多元）同时进行分析时，常常假设因变量组合成的向量服从一个多元正态分布。比如重复测量数据将重复的测量结果（比如各个时间点上的测量结果）视为不同的因变量，可以采用多元方差分析，此时就要求各个因变量的组合向量服从多元正态分布。对多元正态分布的判断通常采用的边际分布来判断，即每个因变量的分布呈正态或近似正态。但实际上单因变量正态是多因变量多元正态的必要非充分条件：所有因变量的组合服从多元正态分布，每个因变量的分布（边际分布）必然呈正态；每个因变量的分布呈正态分布，所有因变量的组合未必呈正态分布。只要有一个因变量不服用正态分布则组合分布肯定不服从多元正态分布。

虽然我们可以通过边际分布来大体判断，有没有统计方法直接判断是否满足多元正态分布呢？当然，R有无所不能的“包”！比如程序包mvnormtest【函数mshapiro.test】、程序包mvShapiroTest【mvShapiro.Test】、程序包MVN【函数mvn】、程序包mvnTest【函数AD.test、CM.test、DH.test、HZ.test、R.test、S2.test等】、程序包mvtnorm【函数pmvnorm】、程序包energy【函数mvnorm.e、mvnorm.test、mvnorm.etest】……

我们以mshapiro.test {mvnormtest}和mvn {MVN}为例演示多元正态分布检验的R实现方法。

mshapiro.test {mvnormtest}

mshapiro.test(U)：Performs the Shapiro-Wilk test for multivariate normality。U要求是数值型矩阵，且a matrix with number of columns (sample size) between 3 and 5000。因为在一般在录入时行表示观测/记录，列表示变量/字段，因此使用此函数需要先对行和列进行转置，转置函数为t(“矩阵或数据框”)。

示例：两组新生儿出生时的体重与身长数据如下。预对两组的体重和身高同时进行检验，可以考虑多元方差分析，应考察多元正态性。

文件导入：
library(openxlsx) #调用程序包openxlsx
mn<-read.xlsx(“D:/Temp/multivnorm.xlsx”,1) #从名称为multivnorm的excel文件中导入第1个sheet的数据到数据框mn中

**将数据调整成mshapiro.test的分析格式：**对数据框mn的行列进行转置
A<-t(mn[1:8,3:4])
B<-t(mn[9:16,3:4])

多元正态分析：
library(mvnormtest) #调用程序包mvnormtest
mshapiro.test(A)
mshapiro.test(B)

结果显示A组多元正态分布SW检验W=0.912，P
=0.366>0.05，呈二元正态分布；B组多元正态分布SW检验W=0.861，P=0.122>0.05，呈二元正态分布。

mvn {MVN}

mvn(data, subset = NULL, mvnTest = c(“mardia”, “hz”, “royston”, “dh”, “energy”), covariance = TRUE, tol = 1e-25, alpha = 0.5, scale = FALSE, desc = TRUE, transform = “none”, R = 1000, univariateTest = c(“SW”, “CVM”, “Lillie”, “SF”, “AD”), univariatePlot = “none”, multivariatePlot = “none”, multivariateOutlierMethod = “none”, bc = FALSE, bcType = “rounded”, showOutliers = FALSE, showNewData = FALSE)。Data为矩阵或数据框。

示例：10名肥胖患者在医生指导下服用药物减肥，按统一标准记录服药前和服药后1-4周的体重。
这是一个典型的无对照的重复测量数据，如将各个时点的体重视为不同的因变量，因变量（结果变量）就不止一个，此时可考虑多元方差分析（MANOVA）对它们同时进行分析。多元方差分析要求多元正态性、组间方差-协方差矩阵同质性、各因变量间有一定的相关性（个体内不独立，但个体间独立）。

library(openxlsx) #调用程序包openxlsx
mn2<-read.xlsx(“D:/Temp/multivnorm.xlsx”,2) #从名称为multivnorm的excel文件中导入第2个sheet的数据到数据框mn2中
U<-mn2[1:10,2:6]
library(MVN)
mvn(U, mvnTest = c(“dh”), multivariatePlot = “qq”)

注：本例10行6列，列名为id、W0、W1、W2、W3、W4。分析时仅筛选列，行不进行筛选，U<-multivnorm[1:10,2:6]可以直接U<-multivnorm[2:6]或者U<-multivnorm[-1]都可以。

结果：本例采用Doornik-Hansen了多变量正态性检验，当然也可以在语句中修改命令，换成Marida, Royston, Henze-Zirkler’s, E-Statistics等方法。结果显示E=6.79，P=0.75>0.05，数据满足多元正态分布。同时结果还给出了每个变量的Shapiro-Wilk检验结果，结果显示W0、W1、W2、W3、W4均满足正态分布。另外本例同时给出了多元正态性的QQ图。

转自个人微信公众号【Memo_Cleon】的统计学习笔记：多元正态分布检验的R实现方法。

生成随机数：

# 生成随机数
set.seed(1230)      # 随机数种子
y1 <- rnorm(100);  # 标准正态分布N(0,1)
y2 <- rexp(100,2);  # 参数为2的指数分布Exp(2)
y3 <- rt(100,1);  # 自由度为1的t分布t(1)
y4 <- -y2;    # -Exp(2)

1 一元正态的评估

1.1 图像法

1.1.1 直方图

直方图是最简单最直白的方法，但是不太适用于少量数据。

# ----------------- 直方图--------------------- #
par(mfrow=c(2,2))   # par是图像参数设置函数，mfrow=c(2,2)是生成2x2的子图
hist(y1,main='Histogram for 100 random numbers from N(0,1)',xlim=c(-5,5))
hist(y2,main='Histogram for 100 random numbers from Exp(2)')
hist(y3[abs(y3) < 5],main='Histogram for 100 random numbers from t(1)',xlim=c(-5,5),xlab='y3')
hist(y4,main='Histogram for 100 random number fram -Exp(2)')

结果：

1.1.2 Q-Q图

Q-Q图的横坐标是理论的分位数，纵坐标是样本分位数，如果样本服从正态分布，那么图中的点应该与直线大致重合。

# ----------------- QQ图--------------------- #
par(mfrow=c(2,2))
qqnorm(y1,main='Q-Q plot for y1 from N(0,1)');qqline(y1)    # 在Q-Q图中添加直线
qqnorm(y2,main='Q-Q plot for y2 from Exp(2)');qqline(y2)
qqnorm(y3,main='Q-Q plot for y3 from t(1)');qqline(y3)
qqnorm(y4,main='Q-Q plot for y4 form -Exp(2)');qqline(y4)

结果：

那么应该如何分析Q-Q图呢？

因此上图从上到下，从左到右依次为：薄尾(但是程度很弱，基本是直线，可以认为是正态)、右偏、厚尾、左偏。

1.2 峰度和偏度

正态分布的偏度应该是0，峰度应该是3.

# ----------------- 偏度和峰度--------------------- #
library(moments)
# skewness
skewness(y1);kurtosis(y1)
skewness(y2);kurtosis(y2)
skewness(y3);kurtosis(y3)
skewness(y4);kurtosis(y4)

结果：

> skewness(y1);kurtosis(y1)
[1] -0.02211832
[1] 2.451047
> skewness(y2);kurtosis(y2)
[1] 1.572117
[1] 4.872045
> skewness(y3);kurtosis(y3)
[1] -1.896557
[1] 12.09301
> skewness(y4);kurtosis(y4)
[1] -1.572117
[1] 4.872045

第一个比较接近标准正态的偏度与峰度，而其他的都明显不符合。

1.3 统计检验

以上方法都存在明显缺陷，图像方法对于小样本并不适用，图像方法以及偏度峰度只提供了一个粗糙而不正式的检验方法，没有一个明确的决定准则。

因此接下来做正式的统计检验，他们基于以下假设：

$H_0:$数据来自正态分布

$H_1:$数据不来自正态分布

1.3.1 Shapiro-Wilks检验

Shapiro-Wilks检验统计量为：

$$W=\frac{{\left({\sum }_{i=1}^n{a}_{(i)}{y}_{(i)}\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\bar{y}\right)}^2}$$

这里$y_{(i)}$是第$i$个样本次序统计量；$a_{(i)}$是标准正态分布中第$i$个次序统计量标准化的期望值。该方法有三个值得注意的地方：

• $\sqrt{W}$约等于实际数据与正态得分之间的相关系数；
• 当$W=1$时，数据恰好完全是正态分布；
• “$W$显著小于1”则表明数据非正态。

代码：

# ----------------- Shapiro-Wilks检验--------------------- #
shapiro.test(y1)  # p越小说明在原假设成立的条件下该事件越不容易发生
shapiro.test(y2)
shapiro.test(y3)
shapiro.test(y4)

结果：

> shapiro.test(y1)  # p越小说明在原假设成立的条件下该事件越不容易发生

	Shapiro-Wilk normality test

data:  y1
W = 0.98635, p-value = 0.3952

> shapiro.test(y2)

	Shapiro-Wilk normality test

data:  y2
W = 0.80837, p-value = 4.513e-10

> shapiro.test(y3)

	Shapiro-Wilk normality test

data:  y3
W = 0.82433, p-value = 1.488e-09

> shapiro.test(y4)

	Shapiro-Wilk normality test

data:  y4
W = 0.80837, p-value = 4.513e-10

以上数据中，y1的$W$值非常接近于1，且$p$值很大，因此接受原假设。而其他几组数据的$p$非常小，而我们知道，$p$越小，可以认为在原假设成立的条件下该样本越不可能发生，在这里就是越不可能是正态分布的数据。

1.3.2 Kolmogorov-Smirnov 检验

Kolmogorov-Smirnov检验的统计量为：

$$D=\sqrt{n}\underset{y}{\sup }\left|F_n(y)-F_0(y)\right|$$

这里的$F_n(y)$是数据的经验累积分布函数(cdf)；$F_0(y)$是与数据同均值、同方差的正态分布的累积分布函数。

该检验方法值得注意的一点是：若$D$值很大，则拒绝原假设$H_0$，也就是$D$越小越好，怎么理解呢？当数据很接近正态分布时，那么按理来说，$F_n(y)$应该比较接近$F_0(y)$，$D$也就越小。

代码：

# ----------------- Kolmogorov-Smirnov检验--------------------- #
library(nortest)    # package for normality test
ks.test(y1,'pnorm',mean(y1),sqrt(var(y1)))
ks.test(y2,'pnorm',mean(y2),sqrt(var(y2)))
ks.test(y3,'pnorm',mean(y3),sqrt(var(y3)))
ks.test(y4,'pnorm',mean(y4),sqrt(var(y4)))

结果：

> ks.test(y1,'pnorm',mean(y1),sqrt(var(y1)))

	One-sample Kolmogorov-Smirnov test

data:  y1
D = 0.043044, p-value = 0.9925
alternative hypothesis: two-sided

> ks.test(y2,'pnorm',mean(y2),sqrt(var(y2)))

	One-sample Kolmogorov-Smirnov test

data:  y2
D = 0.1882, p-value = 0.001676
alternative hypothesis: two-sided

> ks.test(y3,'pnorm',mean(y3),sqrt(var(y3)))

	One-sample Kolmogorov-Smirnov test

data:  y3
D = 0.15339, p-value = 0.01809
alternative hypothesis: two-sided

> ks.test(y4,'pnorm',mean(y4),sqrt(var(y4)))

	One-sample Kolmogorov-Smirnov test

data:  y4
D = 0.1882, p-value = 0.001676
alternative hypothesis: two-sided

第一个的$p$值很大，而其他的都要小于0.05，因此第一个接受原假设，其他的拒绝原假设。

此外，K-S检验还常常用于其他分布的检验，举个例子，为验证每小时参与讨论的人数是否服从泊松分布，原理与正态分布检验是一致的。可以用MATLAB软件进行K-S检验，主要代码程序如下：

lambda = poissfit(A,alpha);   % 对A进行泊松分布拟合
p = poisscdf(A,lambda);   % 生成经验分布函数
H = kstest(A,[A,p],alpha);  % K-S检验

其中A表示由每小时讨论人数组成的向量，alpha为置信率。如果运行结果H=0，则表示A服从泊松分布。对热点话题样本运行的结果显示H=0，因此每小时参与讨论的人数服从泊松分布。

1.3.3 Cramer-von Mises检验

Cramer-von Mises检验的检验统计量为：
$$C=\int {\left[F_n(y)-F_0(y)\right]}^{2}d{F}_0(y)$$

$C$也是和以上的几个检验差不多，越小越好，原理同Kolmogorov-Smirnov检验差不多。

代码：

# ----------------- Cramer-von Mises检验--------------------- #
library(nortest)
cvm.test(y1)
cvm.test(y2)
cvm.test(y3)
cvm.test(y4)

结果：

> cvm.test(y1)

	Cramer-von Mises normality test

data:  y1
W = 0.025708, p-value = 0.8987

> cvm.test(y2)

	Cramer-von Mises normality test

data:  y2
W = 1.1379, p-value = 7.37e-10

Warning message:
In cvm.test(y2) :
  p-value is smaller than 7.37e-10, cannot be computed more accurately
> cvm.test(y3)

	Cramer-von Mises normality test

data:  y3
W = 0.70914, p-value = 5.184e-08

> cvm.test(y4)

	Cramer-von Mises normality test

data:  y4
W = 1.1379, p-value = 7.37e-10

Warning message:
In cvm.test(y4) :
  p-value is smaller than 7.37e-10, cannot be computed more accurately

结果分析同上，第一个$p$很大，因此接受原假设，而其他的则拒绝。

1.3.4 Anderson-Darling检验

Anderson-Darling检验的检验统计量为：

$$A=\int \frac{{\left[F_n(y)-F_0(y)\right]}^2}{F_0(y)\left(1-F_0(y)\right)}d{F}_0(y)$$

对这个检验的理解和上面Kolmogorov-Smirnov检验是一样的，就不多说了。

代码：

# ----------------- Anderson-Darling检验--------------------- #
ad.test(y1)
ad.test(y2)
ad.test(y3)
ad.test(y4)

结果：

> ad.test(y1)

	Anderson-Darling normality test

data:  y1
A = 0.24449, p-value = 0.7564

> ad.test(y2)

	Anderson-Darling normality test

data:  y2
A = 6.5289, p-value = 3.969e-16

> ad.test(y3)

	Anderson-Darling normality test

data:  y3
A = 4.0526, p-value = 3.724e-10

> ad.test(y4)

	Anderson-Darling normality test

data:  y4
A = 6.5289, p-value = 3.969e-16

结果理解同上，不多说了。

2 多元正态分布的评估

有三种方法来检验一个$p$维总体$y$的随机样本 { y 1 , ⋯   , y n } \{y_1,\cdots,y_n \} {y1​,⋯,yn​}是否来自于多元正态分布：

• 1、检验向量的每一维是否都是一元正态分布——如果数据服从多元正态分布，那么每一维都服从一元正态分布，如果检验出有一个不是正态分布，那么联合分布就不是多元正态；
• 2、检验是否每一组二维散点图都没有线性趋势——如果数据服从多元正态分布，那么两个维度之间要么独立，要么成线性关系；
• 3、根据QQ图，检验统计距离 { d 1 , ⋯   , d n } \{d_1,\cdots,d_n\} {d1​,⋯,dn​}是否距离${\chi }^2(p)$很远，其中统计距离定义为：
  $$d_i={\left({\mathbf{y}}_i-\overline{\mathbf{y}}\right)}^{\prime }{\mathbf{S}}^{-1}\left({\mathbf{y}}_i-\overline{\mathbf{y}}\right)$$
  注意，这只是一种近似方法。

以下以一个例子讲解上面三种方法。

2.1 一元检验

首先我们检查每一个变量的QQ图：

二氧化硫分布比较集中，降水以及降水天数背离了正态性；制造企业数和人口数存在很多异常值。综上，该多元数据不服从多元正态分布。

2.2 线性关系检验

接着，绘制两两散点图矩阵：

非线性部分显示了数据与多元正态分布的偏离，显然也说明不是多元正态。

2.3 多元QQ图检验

进一步地，我们绘制整体QQ图：

该图除了检验正态性这一用处外，也可以用来发现可能的异常值。

如果正态性不成立，可以采用一些变量变换方法来获取正态性，如Box-Cox 变换（见我之前的文章，有MATLAB实现方法及原理）。

2.4 R语言实现

# ----------------- 多元正态性--------------------- #
library(HSAUR2)   # which contains the data
# 每一个的QQ图(individual Q-Q plot)
layout(matrix(1:8,4,2,byrow = TRUE))    # 4x2最好不要改动了，不然会出现错误： Error in plot.new() : figure margins too large 
sapply(colnames(USairpollution), function(x)
  {
  qqnorm(USairpollution[[x]], main=x);
  qqline(USairpollution[[x]])
})

# Scatterplot matrix
pairs(USairpollution)

# chi-square Q-Q plot with outlier detection
par(mfrow=c(1,1))
y <- USairpollution
cm <- colMeans(y)
S <- cov(y)
d <- apply(y,1,function(y) t(y-cm) %*% solve(S) %*% (y-cm))
plot(qc <- qchisq( (1:nrow(y)-1/2)/nrow(y),df=7),
     sd <- sort(d),
     xlab = expression(paste(chi[7]^2,"Quantile")),
     ylab = "Ordered distances",
     xlim = range(qc)*c(1,1.1))
oups <- which(rank(abs(qc - sd), ties="random") > nrow(y) - 3)
text(qc[oups],sd[oups]-1.5,names(oups))
abline(a=0,b=1)

对以上代码的理解：

• layout(matrix(c(1,1,2,3),2,2,byrow=TRUE))，其中的2，2表示2×2的图形矩阵；byrow=TRUE表示按列排图，其中的c(1,1,2,3),中有3个不同的值1，2，3，表示2×2的图形矩阵中有3个图，按
  1 1
  2 3
  进行排列。其中1 ，1表示图hist(wt)位于图形矩阵的第一行的1，2列；2表示图hist(mpg)位于2行1列；3表示图hist(disp)位于2行2列。得到上面的图形矩阵。

• sapply()类似于MATLAB中的arrayfun()函数，可以对一个向量批量使用某一个函数。

• pairs()是绘制两两的散点图，pairs在英文中就代表着成双成对，很好理解。

• par()是修改图形参数的函数，与layout()作用差不多。

• apply()函数的使用方法可以见apply函数的使用方法。

• which()函数返回的是满足相应条件的下标/索引，使用方法见which函数的使用方法。

• rank()函数是获得每一个元素排序后的序号，“random” 是相同元素随机编排次序，避免了“先到先得”，“权重”优于“先后顺序”的机制增大了随机的程度。使用方法见rank函数的使用方法。

• text()函数就是添加文本标记，它是向绘图区域内部添加文本，而mtext()则向图形的四个边界之一添加文本。

• abline()则常用用于为图形添加参考线，其使用格式为abline(h=yvalues, v=xvalues, params)，其中h为水平方向添加，而v为纵向。

练习题

练习1

代码：

# 练习
data <- c(-0.6,3.1,25.3,-16.8,-7.1,-6.2,25.2,22.6,26.0)
# (a)
qqnorm(data,mian='1996-2005年间道琼斯工业平均指数的年回报率');qqline(data)

# (b)
shapiro.test(data)

结果：

> shapiro.test(data)

	Shapiro-Wilk normality test

data:  data
W = 0.85607, p-value = 0.08688

练习2

代码：

Y1 <- c(1,2,3,3,4,5,6,8,9,11)
Y2 <- c(18.95,19.00,17.95,15.54,14.00,12.95,8.94,7.49,6.00,3.99)
# (a)
data <- data.frame(Y1,Y2)
mu <- c(0,0)
mu[1] <- mean(Y1)
mu[2] <- mean(Y2)
cov.y1y2 <- cov(data)
cor.y1y2 <- cor(data)

# (b)
d <- apply(data,1,function(data) t(data-mu) %*% solve(cov.y1y2) %*% (data-mu))

# (c)
plot(qc <- qchisq((1:length(d) - 1/2)/length(d) , df=length(d)-1),
     sd <- sort(d))

结果：

# (a)
> mu
[1]  5.200 12.481
> cov.y1y2
          Y1        Y2
Y1  10.62222 -17.71022
Y2 -17.71022  30.85437
> cor.y1y2
           Y1         Y2
Y1  1.0000000 -0.9782684
Y2 -0.9782684  1.0000000

# (b)
> d
 [1] 1.8753045 2.0203262 2.9009088 0.7352659 0.3105192 0.0176162 3.7329012 0.8165401 1.3753379 4.2152799

( C ):

附录——对p值的理解

求法方面的理解:p值可以理解为在原假设成立的情况下当前样本或倾向于备择假设成立的极端情况的发生概率

含义方面的理解:在原假设成立的情况下拒绝原假设的最小显著性水平。如果显著性水平α>p，拒绝原假设，反过来则接受原假设

九、独立性检验和正态性检验

1.独立性检验

独立性检验，指的是将一个多元总体$X\sim N_p(\mu ,\Sigma )$划分成$k$个部分，探究每个部分之间是否独立的问题，这样做的好处是显而易见的，如果一个总体$X$可以划分成多个独立的部分，那么只需要对每一个部分分开讨论即可，无疑降低了运算量。在多元统计中，可以视为有如下分解：
$$X=\left[\begin{array}{c}{X}^{(1)}\\ ⋮\\ X^{(k)}\end{array}\right],\mu =\left[\begin{array}{c}{\mu }^{(1)}\\ ⋮\\ {\mu }^{(k)}\end{array}\right],\Sigma =\left[\begin{array}{ccc}{\Sigma }_{11}& \cdots & {\Sigma }_{1k}\\ ⋮& & ⋮\\ {\Sigma }_{k1}& \cdots & {\Sigma }_{kk}\end{array}\right].$$
每一个分向量$X^{(t)}$都是$p_t$维的，对应的${\mu }^{(t)}$也是$p_t$维的，${\Sigma }_{tt}$是$p_t\times p_t$的。在多元正态分布的介绍中提到，如果$X^{(1)},\cdots ,X^{(k)}$是独立的，那么${\Sigma }_{ij}=O$对任何$i\ne j$都成立，反之也成立，因此在正态总体下，假设检验就变成了以下的形式：
$$H_0:\forall i\ne j,{\Sigma }_{ij}=O\iff H_1:\exists i\ne j,{\Sigma }_{ij}\ne O.$$
由于样本均值、样本离差阵是对总体均值、自协方差矩阵的估计，因此我们也可以对样本均值和样本离差阵作同型分解。如果$H_0$成立，则$X_{(\alpha )}^{(t)}\sim N_p({\mu }^{(t)},{\Sigma }_{tt})$且相互独立，那么似然函数就是
$$L(\mu ,\Sigma )=\prod _{t=1}^k{L}_t({\mu }^{(t)},{\Sigma }_{tt}),$$
取最大值的情况显然是${\mu }^{(t)}={\bar{X}}^{(t)},{\Sigma }_{tt}=A_{tt}/n$，所以似然比统计量的分子是
$$\begin{array}{rl}& \prod _{t=1}^n(2\pi {)}^{-n{p}_t/2}|A_{tt}/n{|}^{-n/2}\exp \left\{-\frac{1}{2}\sum _{\alpha =1}^n(X_{(\alpha )}^{(t)}-{\bar{X}}^{(t)}{)}^{\prime }{\left(\frac{A_{tt}}{n}\right)}^{-1}(X_{(\alpha )}^{(t)}-{\bar{X}}^{(t)})\right\}\\ =& (2\pi {)}^{-np/2}\exp \left\{-\frac{1}{2}\sum _{\alpha =1}^n(X_{(\alpha )}-\bar{X}{)}^{\prime }{\left(\frac{A}{n}\right)}^{-1}(X_{(\alpha )}-\bar{X})\right\}\prod _{t=1}^k{\left|\frac{A_{tt}}{n}\right|}^{-n/2}.\end{array}$$
这里的转换可以用之前常用的迹变换得出。观察分子与分母，发现其大部分是相同的，所以得到似然比统计量为
$$\lambda =\frac{{\prod }_{t=1}^k|A_{tt}/n{|}^{-n/2}}{|A/n{|}^{-n/2}}={\left(\frac{|A|}{{\prod }_{t=1}^n|A_{tt}|}\right)}^{n/2}\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{V}^{n/2}.$$
所以我们取检验统计量为
$$V=\frac{|A|}{{\prod }_{i=1}^k|A_{tt}|}.$$
并且有结论保证，在$H_0$成立的条件下，$-b\ln V\stackrel{H_0}{\to }{\chi }^2(f)$，这里
$$\begin{gathered} b=n-\frac{3}{2}-\frac{p^3-{\sum }_{t=1}^k{p}_t^3}{3(p^2-{\sum }_{t=1}^k{p}_t^2)}, \\ f=\frac{1}{2}\left[p(p+1)-\sum _{t=1}^k{p}_t(p_t+1)\right]. \end{gathered}$$
事实上$-b\ln V$是$-2\ln \lambda$的近似，故$b$也是$n$的近似，而$f$就是两个参数空间的维度之差。

2.一元数据正态性检验

回顾我们之前提到的假设检验，包括均值向量、自协方差矩阵、独立性的检验，都基于一个前提——总体是多维正态分布，如果这个正态性不满足，与三大分布相关的统计量转化、似然比统计量的表现形式都将不同于此形式，从而无法应用已有的结论。因此，本节探讨样本的正态性检验，概括起来就是，给定$n$个$p$维样本$X_{(\alpha )}$，判断总体$X$是否服从$N_p(\mu ,\Sigma )$分布。

多元数据的正态性检验问题，常常转化为多个一元或二元数据的正态性检验，或者先求$X$的分量的线性组合再化为一元数据的正态性检验等。虽然我们知道，边缘分布的正态性不能推出总体分布的正态性，但是在实际应用中，这种情况并不常见，所以我们可以先将目光放在一元数据的正态性检验。

常用于一元数据检验的方法有Pearson ${\chi }^2$检验法（比较适合离散情形）、Kolmogorov检验法（比较适合连续情形），不过在Kolmogorov检验中我们需要得知总体的参数，即均值和方差，在实际应用中这个条件很难满足，所以我们会使用总体均值和总体方差代替，这就是Lilliefors检验。

还有一些仅适用于正态分布的检验法：偏度峰度检验法，Q-Q图和P-P图检验法、Anderson-Darling统计量检验法、Cramer-von Mises统计量检验法等。

偏度峰度法指的是，计算样本偏度和样本峰度：
$$G_1=\frac{\sum (X_i-\bar{X}{)}^3}{[\sum (X_i-\bar{X}{)}^2{]}^{3/2}},G_2=\frac{\sum (X_i-\bar{X}{)}^4}{[\sum (X_i-\bar{X}{)}^2{]}^2},$$
在正态性成立时，近似有
$$\begin{gathered} G_1\sim N\left(0,\frac{6(n-2)}{(n+1)(n+3)}\right), \\ {G}_2\sim N\left(3-\frac{6}{n-1},\frac{24n(n-2)(n-3)}{(n+1{)}^2(n+3)(n+5)}\right). \end{gathered}$$
很容易用Z检验找到其拒绝域。

Q-Q(Quantile Quantile)图检验法是一种图示检验法，绘制$(q_i,x_{(i)}^{\ast })$散点图，这里$q_i={\Phi }^{-1}(p_i)$是样本的$p_i$分位数，$x_{(i)}^{\ast }$是样本的$p_i$分位数，如果$X$是一元正态总体，则这些散点应该散布在一条直线上。P-P图检验法也是图示检验，绘制的数据点是$(p_i,F(x_{(i)}^{\ast }))$，其中$p_i$是经验分布函数$F_n(x)$在$x_{(i)}^{\ast }$上的值，$F(x_{(i)}^{\ast })$是$\Phi (x)$在$x_{(i)}^{\ast }$上的值。在实际应用Q-Q图检验和P-P图检验时，$x_{(i)}^{\ast }$要先选好。

Anderson-Darling $A^2$检验（AD检验）的检验统计量是
$$A^2=n{\int }_{-\infty }^{\infty }\frac{(F_n(x)-\Phi (x){)}^2}{\Phi (x)(1-\Phi (x))}\mathrm{d}\Phi (x),$$
这里$[\Phi (x)(1-\Phi (x)){]}^{-1}$是权重函数，如果权重函数取$1$，就得到Cramer-von Mises $W^2$检验的检验统计量
$$W^2=n{\int }_{-\infty }^{\infty }(F_n(x)-\Phi (x){)}^2\mathrm{d}\Phi (x).$$
结合Kolmogorov-Smirnov统计量$D=\sup |F_n(x)-\Phi (x)|$，这三个统计量都是原假设成立时不能过大的，依赖于一个概率表值来检验原假设是否应该被接受。不过，这三种检验方式适用于各种假设检验，只要将表达式中的$\Phi (x)$换成对应的分布函数即可。

3.多元数据的正态性检验

对于二元数据，存在一种粗糙的检验方法：等概椭圆检验法。其理论基础是二维随机向量$X$如果来自于正态总体，则其概率密度函数等高线应该是一个椭圆，即$X\sim N_2(\mu ,\Sigma )$时，应有
$$f(x_1,x_2)=a\iff (X-\mu {)}^{\prime }{\Sigma }^{-1}(X-\mu )=b^2.$$
所以我们计算二元数据$X_{(i)}$到$\bar{X}$的马氏距离$D_i=(X_{(i)}-\bar{X}{)}^{\prime }{S}^{-1}(X_{(i)}-\bar{X})$，在给定数值$p_0$下，$D_i\le p_0$的频率应该和某一个定值比较接近，这个定值可以通过查表获得。由于这是一种比较粗糙的方法，我们在实际应用中会使用更为正式的方法。

现在介绍$p$维数据${\chi }^2$统计量的Q-Q图检验法，我们将假设确定为参数已知的，即
$$H_0:X\sim N_p(\mu ,\Sigma )\iff H_1:X\nsim N_p(\mu ,\Sigma ).$$
由于在正态性假设$H_0$成立的前提下，样本$X$到中心$\mu$的马氏距离存在以下关系：
$$D^2=(X-\mu {)}^{\prime }{\Sigma }^{-1}(X-\mu )\sim {\chi }^2(p),$$
所以我们可以直观地想到验证样本的马氏距离是否具有这样的关系。因此，我们计算样本$X_{(\alpha )}$到$\mu$的马氏距离$D_{\alpha }^2=(X_{(\alpha )}-\mu {)}^{\prime }{\Sigma }^{-1}(X_{(\alpha )}-\mu )$，并对$D_{\alpha }^2$进行排序得到次序统计量$D_{(\alpha )}$，计算其经验分布函数，这样有了经验分布函数与${\chi }^2(p)$分布的分布函数后，就可以绘制Q-Q图或者P-P图。

在实际应用中，我们往往不知道$\mu ,\Sigma$的值，所以会用样本均值$\bar{X}$和样本协方差阵$A/(n-1)$代替，得到的Q-Q图或P-P图应该是一条通过原点、斜率为1的直线，如果是这样，就可以接受正态性假设，否则应当拒绝。

回顾总结

1. 正态总体的独立性检验，我们一般会取检验统计量为
   $$V=\frac{|A|}{{\prod }_{t=1}^k|A_{tt}|}.$$
   当$n\to \infty$时，有$-b\ln V\to {\chi }^2(f)$，这里
   $$\begin{gathered} b=n-\frac{3}{2}-\frac{p^3-{\sum }_{t=1}^k{p}_t^3}{3(p^2-{\sum }_{t=1}^k{p}_t^2)}, \\ f=\frac{p(p+1)}{2}-\sum _{t=1}^k\frac{p_k(p_k+1)}{2}. \end{gathered}$$

2. 一元总体的正态性检验有很多方法，如K-S检验、A-D检验、Cramer-von Mises检验，但K-S检验的效果一般，A-D检验的效果比较好，其检验统计量是
   $$A^2=n{\int }_{-\infty }^{\infty }\frac{(F_n(x)-\Phi (x){)}^2}{\Phi (x)(1-\Phi (x))}\mathrm{d}\Phi (x).$$

3. Q-Q图是分位数图，首先选定一组分位数间隙$x_{(i)}^{\ast }$，然后在样本中寻找相应分位数，在总体中也寻找相应分位数，将分位数绘制成散点图，观察其是否位于一条直线上。

4. P-P图是累计分布图，首先选定一组分位数间隙$x_{(i)}^{\ast }$，然后绘制经验分布函数与总体分布函数在$x_{(i)}^{\ast }$处的取值，将两个取值绘制成散点图，观察其是否位于一条直线上。

5. 多元总体的正态性检验采用${\chi }^2$统计量的Q-Q图检验法，计算样本到中心$\bar{X}$的马氏距离并排序，用Q-Q图判断是否属于${\chi }^2(p)$分布，或用K-M检验法。马氏距离的定义如下：
   $$D_{\alpha }=(X_{(\alpha )}-\bar{X}{)}^{\prime }{S}^{-1}(X_{(\alpha )}-\bar{X}).$$