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，拒绝原假设，反过来则接受原假设