正态分布笔记

分布

表示分布最常用的方法是直方图(histogram)，这种图用于展示各个值出现的频数或概率。频数指的是数据集中的一个值出现的次数。概率就是频数除以样本数量n。

用表示概率的直方图称为概率质量函数(Probability Mass Function, PMF)。这个函数是值到其概率的映射。

连续型变量的分布

连续型变量是在连续区间取值，例如质量、长度、面积、体积、寿命、距离等就是连续型变量，现在试想一下连续变量观测值的直方图；如想其纵坐标为相对频数，那么所有这些矩形条的高度和为1，那么完全可以重新设置量纲，例这些矩形条的面积为1，如果不断增加观测值，并不断增加直方图的矩形条的数目，这些直方图就会越来越像一条光滑曲线，其下面的面积和为1，这种曲线就是概率密度函数(probability density function, pdf)，简称为密度函数或密度，下图就展示了逐渐增加矩形条数的直方图和一个形状类型的密度曲线：

par(mfrow=c(2,2))
x=rnorm(10000)
z<-seq(-4,4,length=10000)
hist(x,20,probability=T,col=2)
hist(x,50,probability=T,col=2)
hist(x,100,probability=T,col=2)
plot(z,dnorm(z),type="l",
     xlab="x",ylab="Density Function",
     bty="n",
     main="Dencity Function of x")

什么是密度？

这里说一下密度，密度是描述变量分布的一个重要概率，简单地说，它描述变量取各个值的可能性大小，变量在密度曲线高的位置取值的可能性比在低的位置取值的可能性大。在概率密度曲线中，纵坐标就是密度(Density)。

而密度曲线(density curve)是分布的一个数字模型，按照定义，任何密度曲线下的总面积都是1，在一个区间的密度曲线下的面积为出现在这个区间中的观测值占所有观测值的比例。

在矩形面积和为1的直方图中，连续变量落入某个矩阵的概率为该矩形的面积，类似地，连续变量落入一个区间的概率就是概率密度函数曲线在这个区间上所覆盖的面积，也就是密度函数在这个区间上的积分，根据初等微积分，连续函数在一个点的积分是0（因为曲线下面的面积退化成一条垂直于x轴，平行于y轴的直线），因此对于连续变量，取任何一个单独特定值的概率都是0，而只有变量取值于某个（或若干个）区间的概率才有可能大于0。连续变量和离散变量类似，它的密度函数（这时用f表示随机变量X的密度函数）必须满足$f(x) \geq 0$，而且密度曲线下面覆盖的总面积为1，即：

$\int_{-\infty} ^{ \infty }{f(x)}\,{\rm d}x=1$

而X落入区间(a,b)的概率则为：
$P(a\leq X \leq b)=\int_a^b{f(x)\,{\rm d}x}$

连续变量也有累积分布函数(cumulative distribution function, CDF)，定义为：

$F(x)=P(X\leq x)=\int_{-\infty} ^{x }{f(x)}\,{\rm d}x$

这样，和离散变量一样，如果要求P(a<X≦b)的值，就可以通过两个累积分布F(b)和F(a)之间在用公式获得，如下所示：

$P(a\leq X \leq b)=P(X\leq b)-P(X\leq a)=F(b)-F(a)$

而不用根据密度函数求积分了。

连续随机变量也有描述变量“位置”的总体均数、总体中位数、总体百分位数以及描述变量分布（集中）程度的总体标准差和总体方差等概念，对于连续分布随机变量X，它的总体均值（通常用μ表示）定义为：
$\mu=E(X)=\int_{-\infty} ^{ \infty }{xf(x)}\,{\rm d}x$

这也称为X的数学期望，变量X的总体方差（通常用$\delta^2$表示），定义为：

$\delta^2=Var(X)=\int_{-\infty} ^{ \infty }{(x-\mu)f(x)}\,{\rm d}x$

总体标准差定义为总体方差的平方根，如下所示：

$\delta=\sqrt{Var(X)}=\sqrt{\int_{-\infty} ^{ \infty }{(x-\mu)f(x)}\,{\rm d}x}$

正态分布

下图是四个二项分布图，这四个二项分布有共同的p=0.1，但是有不同的n，分别为5，20，50，500，可以看出，当n较小时，分布很不对称，当n增加时，分布变得对称，实际上，当n充分大时，二项分布趋向于一个连续型对称分布，这就是正态分布(normal distribution)，如下所示：

n<-c(5,20,50,500)
p<-0.1
par(mfrow=c(2,2))
barplot(dbinom(0:n[1],n[1],p),main="Binomial(5,0.1)")
barplot(dbinom(0:(2*floor(2*n[2]/5)),n[2],p),main="Binomial(20,0.1)")
barplot(dbinom(0:(floor(n[3]/3)),n[3],p),main="Binomial(50,0.1)")
barplot(dbinom(25:75,n[4],p),main="Binomial(500,0.1)")

在实际生活中，近似地服从正态分布的变量也很常见，例如测量误差，商品的重量或尺寸，某特定人群的身高和体重等，此外，在一定条件下，许多不是正态分布的样本均值在样本量很大时，也可用正态分布来近似估计。

正态分布也是一族分布，不同的正态分布根据它们的均值和标准差不同而有所区别，因此一个正态分布用$N(\mu,\delta)$表示，其中$\mu$是总体均值，而$\delta$是总体标准差，正态分布也常用$N(\nu,\delta^2)$来表示，这里的$\delta^2$为总体方差，这些总体参数在实际问题中是未知的，但可以估计，例如用样本均值和样本标准差来估计总体均值和总体标准差，正态分布的密度曲线是一个对称的钟形曲线（最高点在均值处），正态分布$N(\mu,\delta)$的密度函数如下所示：

正态分布的推导

在理解正态分布的时候，如果知道它的推导过程，以及曲线下面积是如何求的，就会对正态分布的理解更深入，但是这涉及到一定的数学知识，包括大数定量，微积分等，目前我没有能力理解这些，具体的可以看陈希孺的《概率论与数理统计》等资料。不过真正用这个函数来计算正态变量概率数值的情况不多，多数直接采用累积分布函数来计算的。

正态分布密度曲线

下面分别绘制均值（以下用m表示）为0，标准差（以下用sd表示）为1，m=0，sd=0.5，m=0，sd=2，m=-2，sd=1的正态分布密度曲线（密度曲线下面的积分面积对应的就是相应范围内的概率），可以看出不同的均值与标准差，正态分布的密度曲线形状就不一样：

set.seed(1000)
x <- seq(-5,5,length.out=100)
y <- dnorm(x,0,1)
plot(x,y,col="red",xlim=c(-5,5),ylim=c(0,1),type="l",
xaxs="i",yaxs="i",ylab="density",xlab="",
main="The Normal Density Distribution")
lines(x,dnorm(x,0,0.5),col="green")
lines(x,dnorm(x,0,2),col="blue")
lines(x,dnorm(x,-2,1),col="orange")
legend("topright",legend=paste("m=",c(0,0,0,-2),"sd=",c(1,0.5,2,1)),lwd=1,col=c("red","green","blue","orange"))

结果如下所示：

正态分布的概率计算

作为一个连续变量，正态变量落在某区间的概率不是等于在这个区间上，密度曲线下面的面积，下图的阴影面积等于标准正态分布变量落在区间(-1.5,1.5)中的概率，如果密度函数为$\phi(x)$，那么这个面积就等于以下的积分，如下所示：
$\int_{-1.5}^{1.5}{\phi(x)}\,{\rm d}x=\Phi(1.5)-\Phi(-1.5)$

在计算机未广泛应用之前，为了方便计算，数学家们就编制了标准正态分布表，而在计算机广泛使用之后，则很容易计算出服从正态分布的数据的某个范围的概率，但为了理解正态分布的原理，还是先看一下如何使用标准正态分布表来计算一下标准正态分布变量落在区间(-1.5,1.5)中的概率。

在使用标准正态表之前，我们知道这些信息：

1. 标准正态分布是关于y轴对称的；
2. 标准正态分布的概率是落在-1.5到1.5这个区间内的曲线下面积；

现在查表，可以知道1.5对应的这个数值是0.9332，它表示小于1.5的曲线下面积是0.9332，-1.5没法查，但是可以查-1.5对称的曲线下面积，-1.5关于y轴对称的直线是0.5，它的曲线下面积为0.6915，也就是说小于0.5的曲线下面积为0.6915，那么大于0.5的曲线下面积为1-0.6915=0.3085，这个0.3085其实就是小于-1.5的曲线下面积，那么标准正态分布变量落在区间(-1.5,1.5)中的概率就是0.9332-0.3085=0.6247

也就是这下面的这个灰色面积，如下所示：

x<-c(seq(-4,4,length=1000))
y<-dnorm(x)
r1<-0.5
r2<-1.5
x2<-c(r1,r1,x[x<r2&x>r1],r2,r2)
y2<-c(0,dnorm(c(r1,x[x<r2&x>r1],r2)),0)
plot(x,y,type="l",main="P(-0.5<Z<1.5)=0.6246553",xlab="z",ylab=expression(phi(z)))
abline(0,0)
polygon(x2,y2,col="grey")

不过，现在一般不用查表，使用软件计算即可，如下所示：

> pnorm(1.5)-pnorm(-0.5)
[1] 0.6246553

正态分布累积曲线

代码如下所示：

set.seed(1000)
x <- seq(-5,5,length.out=100)
y <- pnorm(x,0,1)
plot(x,y,col="red",xlim=c(-5,5),ylim=c(0,1),type="l",
xaxs="i",yaxs="i",ylab="density",xlab="",
main="The Normal Cumulative Distribution")
lines(x,pnorm(x,0,0.5),col="green")
lines(x,pnorm(x,0,2),col="blue")
lines(x,pnorm(x,-2,1),col="orange")
legend("bottomright",legend=paste("m=",c(0,0,0,-2),"sd=",c(1,0.5,2,1)),lwd=1,col=c("red","green","blue","orange"))

结果如下所示：

总体分位数和尾概率

对于连续型随机变量X，α下侧分位数（又称α分位数，α-quantile）定义为满足关系$P(X\leq x)=\alpha$的数x，这里的α称为下（左）侧尾概率(lower/left tail probability)。

α上侧分位数（又称α上分位数，α-upper quantile）定义为它满足$P(X\geq x)=\alpha$的数x，这里的α称为上（右）侧尾概率(upper/right tail probability)。

对于连续分布，α上侧分位数等于(1-α)下侧分位数，而(1-α)下侧分位数等于α上侧分位数，通常用$z_{\alpha}$表示标准正态分布α的上侧分位数，即对于标准正态分布变量Z，有$P(Z>z_{\alpha})=\alpha$，有些文献用符号$z_{1-\alpha}$而不是$z_{\alpha}$来表示α上侧分位数。

下图表示的就是标准正态分布0.05上侧分位数$z_{\alpha}=z_{0.05}(=1.645)$，图形右边在区间$(z_{0.05},+\infty)$曲线下的面积等于相应的尾概率(α=0.05)：

x<-c(seq(-4,4,length=1000))
y<-dnorm(x)
plot(x,y,type="l",xlab="z",ylab=expression(phi(z)))
abline(0,0)
r<-qnorm(0.05,low=F)
polygon(c(r,r,x[x>r]),c(0,dnorm(c(r,x[x>r]))),col="red")
text(0,0.2,expression(P(Z<z[0.05])==0.95))
text(2.5,0.08,expression(P(Z>z[0.05])==0.05))
text(1.5,0.02,expression(P(z[0.05])==1.645))

随后要介绍的$\chi^2$分布，t分布和F分布都是由正态分布变量导出的三种分布。

正态分布案例之一

再看一个案例：

B县的全体高一学生参加某补习班的数学测验，计算分数后得知“数学测验结果”可视为平均值为45，标准差为10的正态分布，那么这些数据就可以表示为这个样子：

绘制出平均值为45，标准差为10的正态分布表示，下图斜线部分的面积是0.5，如下所示：

还可以得到出这些信息：

第一：得分在45分以上的考生比例，占全体考生总数的0.5（50%）；

第二：从全体考生中，随机抽出一个人，其得分在45分以上的概率为0.5（50%）；

第三，将这些数学成绩标准化后（也就是常说的Z转换），0以上的考比例，占全体考生总数的0.5（50%），如下所示（看红线，不要看0.29）：

再继续这个案例，由于它是符合平均数为45，标准差为10的正态分布，那么，如果我要求某一个分数以上的考生比例，例如63分以上的考生比例，它的计算过程就是先把这个数字转换到标准正态分布上，过程如下所示（这就是Z转换）：

然后从标准正态分布表查1.8，它对应的数字是0.4641，这个数字表示的是小于63这个数字的概率，那么大于63的概率就是：0.5-0.4641=0.0359（3.59%）。它对应的面积就是下面的这个样子：

将它转换为标准正态分布，则是这个样子：

不过现在过计算正态分布并不需要查表，查表只是了解一下这个计算原理，现在都用计算机计算了。

R语言正态分布函数

R语言中与正态分布有关的函数如下所示：

1. 密度函数： dnorm(x, mean = 0, sd = 1, log = FALSE)
2. 分布函数： pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
3. 分位数函数：qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
4. 随机函数： rnorm(n, mean = 0, sd = 1)
   其余的分布函数与正态分布的前缀相同。

正态分布相应的概率计算

求正态N(3.1.5)变量X在区间(2,4)中的概率p，如下所示：

pnorm(4,3,1.5)-pnorm(2,3,1.5)

pnorm(4,3,1.5)-pnorm(2,3,1.5)
[1] 0.4950149

正态分布的检验

判断数据是否符合正态分布通常用的是非参数检验方法（因为总体分布未知）

shapiro.test()函数

Shapiro-Wilk正态分布检验: 用来检验是否数据符合正态分布，类似于线性回归的方法一样，是检验其于回归曲线的残差。该方法推荐在样本量很小的时候使用，样本在3到5000之间。该检验原假设为H0:数据集符合正态分布，统计量W 最大值是1，越接近1，表示样本与正态分布匹配。如果p-value小于显著性水平α(0.05)，则拒绝H0，表示样本不符合正态分布，如果p-value在大于显著性水平，则接受H0，表示样本符合正态分布。

set.seed(1)
s <- rnorm(1000)
shapiro.test(s)

其结果为：
Shapiro-Wilk normality test

data: s
W = 0.9988, p-value = 0.7256

结论: W接近1，p-value>0.05，不能拒绝原假设，所以数据集s符合正态分布。

Kolmogorov-Smirnov连续分布检验

Kolmogorov-Smirnov连续分布检验:检验单一样本是不是服从某一预先假设的特定分布的方法。以样本数据的累计频数分布与特定理论分布比较，若两者间的差距很小，则推论该样本取自某特定分布族。
该检验原假设为H0:数据集符合正态分布，H1:样本所来自的总体分布不符合正态分布。令F0(x)表示预先假设的理论分布，Fn(x)表示随机样本的累计概率(频率)函数.
统计量D为: D=max|F0(x) - Fn(x)|，D值越小，越接近0，表示样本数据越接近正态分布，如果p-value小于显著性水平α(0.05)，则拒绝H0。

Kolmogorov-Smirnov 检验,用法如下： ks.test(x, y, …, alternative = c(“two.sided”, “less”, “greater”), exact = NULL)
R语言中ks.test有四个参数，第一个参数x为观测值向量，第二个参数y为第二观测值向量或者累计分布函数或者一个真正的累积分布函数如pnorm，只对连续CDF有效。第三个参数为指明是单侧检验还是双侧检验，exact参数为NULL或者一个逻辑值，表明是否需要计算精确的P值。

set.seed(1)
s <- rnorm(1000)
ks.test(s,"pnorm")

结果如下：
One-sample Kolmogorov-Smirnov test

data: s
D = 0.0211, p-value = 0.7673
alternative hypothesis: two-sided

D值很小, p-value>0.05，不能拒绝原假设，所以数据集s符合正态分布。

KS的第二案例应用

比较两组数据的分布是否相同

ks.test(rnorm(100),rnorm(50))

Two-sample Kolmogorov-Smirnov test

data: rnorm(100) and rnorm(50)
D = 0.14, p-value = 0.5185
alternative hypothesis: two-sided

在这一案例中，我们比较了两个均值和方差一样的观测值，他们D值很小，p值大于0.05，所以我们不能拒绝两个观测值分布相同的假设。

KS比较某一数据符合某一分布

检验一组数据是否符合均匀分布（生成均匀分布的函数为runif))，如下所示：

ks.test(rnorm(100),"punif")

结果中p-value为2.2e-16，我们用rnorm(100)生成的是正态分布的数据，因此它不符合均匀分布（参数为punif）。如果参数设为pnorm，再看一下。

ks.test(rnorm(100),"pnorm")

结果的p值为0.6722，符合正态分布。

Pearson 拟合优度chi-square检验

（略）

QQ图

正态概率图用于检查一组数据是否服从正态分布。是实数与正态分布数据之间函数关系的散点图。如果这组实数服从正态分布，正态概率图将是一条直线。通常，概率图也可以用于确定一组数据是否服从任一已知分布，如二项分布或泊松分布。

QQ图是quantile—quantile的简写，一个根据正态分布的分位数，一个是样本分位数；如果两者在一条直线上，说明样本符合状态分布的分位数，即服从正态分布。QQ图是一种散点图,对应于正态分布的QQ图,就是由标准正态分布的分位数为横坐标,样本值为纵坐标的散点图. 要利用QQ图鉴别样本数据是否近似于正态分布,只需看QQ图上的点是否近似地在一条直线附近,而且该直线的斜率为标准差,截距为均值. 用QQ图还可获得样本偏度和峰度的粗略信息。

另外，还有一种PP图。P-P图是根据变量的累积比例与指定分布的累积比例之间的关系所绘制的图形。通过P-P图可以检验数据是否符合指定的分布。当数据符合指定分布时，P-P图中各点近似呈一条直线。如果P-P图中各点不呈直线，但有一定规律，可以对变量数据进行转换，使转换后的数据更接近指定分布。 Q-Q图同样可以用于检验数据的分布，所不同的是，Q-Q图是用变量数据分布的分位数与所指定分布的分位数之间的关系曲线来进行检验的。

set.seed(1)
s <- rnorm(1000)
qqnorm(s)
qqline(s,col="red")

结果如下所示：

解释：纵轴是数据的实际分位数，横轴是正态分布的理论分位数，红线是45度的参考线，如果数据是正态分布，则散点应该落在直线上。即实际分位数与理论分位数偏差不大。说明来自正态总体。

峰度与偏度

峰度（Kurtosis）

峰度是描述某变量所有取值分布形态陡缓程度的统计量，它是和正态分布相比较的。
公式如下所示：

$$Kurtosis=\frac{\sum(X-\mu)^{4}}{\delta^{4}} -3$$

正态分布的Kurtosis为K=3,为了描述方便，有时候需要减3（上述公式就减了3），有的公式没有，但原理是相同的。

Kurtosis=0 与正态分布的陡缓程度相同。
Kurtosis>0 比正态分布的高峰更加陡峭——尖顶峰
Kurtosis<0 比正态分布的高峰来得平台——平顶峰

偏度（Skewness）

偏度是描述某变量取值分布对称性的统计量。公式如下

$$Skewness=\sum[(\frac{X-\mu}{\delta})^3]$$

Skewness=0 分布形态与正态分布偏度相同
Skewness>0 正偏差数值较大，为正偏或右偏。长尾巴拖在右边。
Skewness<0 负偏差数值较大，为负偏或左偏。长尾巴拖在左边。

正态分布检验补充案例

words1 = rnorm(100)
words2 = rt(100, df=3)

plot(density(words1))
lines(density(words2))
shapiro.test(words1)
shapiro.test(words2)

par(mfrow=c(2,1))
qqnorm(words1);qqline(words1, col = 2)
qqnorm(words2);qqline(words2, col = 2)

words1与words2密度图如下所示：

从图上大致可以看出来words1与标准正态图比较接近。

用shapiro检验后的结果如下所示：

从计算结果来看，words1是p值为0.09861，大于0.05，不能否定原假定，因此words1符合正态分布。同理,words2不符合正态分布。

QQ图的结果如下所示：

从QQ图中可以看出，words1符合正态分布，而words2不符合正态分布。

R语言中进行正态分布检验的其它包

nortest包

lillie.test()可以实行更精确的Kolmogorov-Smirnov检验。

ad.test()进行Anderson-Darling正态性检验。

cvm.test()进行Cramer-von Mises正态性检验。

pearson.test()进行Pearson卡方正态性检验。

sf.test()进行Shapiro-Francia正态性检验。

fBasics包

ksnormTest() 进行Kolmogorov-Smirnov正态性检验。

shapiroTest() 进行Shapiro-Wilk正态性检验。

jarqueberaTest() 进行jarque-Bera正态性检验。

dagoTest() 进行D’Agostino正态性检验。

adTest() 进行Anderson-Darling正态性检验。

cvmTest() 进行Cramer-von Mises正态性检验。

lillieTest() 进行Lilliefors(KS)正态性检验。

pchiTest() 进行Pearson chi-square正态性检验。

sfTest() 进行Shapiro-Francia正态性检验。

jbTest() Finite sample adjusted JB LM and ALM test。

normalTest() 进行Kolmogorov-Smirnov正态性检验。

Jarque-Bera Test检验正态分布

Jarque-Bera检验是用于总体分布的正态性检验的一种方式。该检验基于数据样本的偏度和峰度，评价给定数据服从未知均值和方差正态分布的假设是否成立，正态分布的偏度（三阶矩）S=0，峰度（四阶矩）K=3，若样本来自正态总体，则他们分别在0，3附近。基于此构造一个包含x2(卡方)统计量：n为样本容量,k为自由度。代码如下所示：

install.packages("tseries")
library(tseries)
set.seed(1)
x <- rnorm(150)
res <- jarque.bera.test(x)
names(res)
res$p.value

运行结果如下所示：

> set.seed(1)
> x <- rnorm(150)
> res <- jarque.bera.test(x)
> names(res)
[1] "statistic" "parameter" "p.value"   "method"    "data.name"
> res$p.value
[1] 0.8601533

p值大于0.05，可以视为符合正态分布。

案例第二：用Jarque-Bera Test检验正态分布

install.packages("Ecdat")
library(Ecdat)
data(SP500)
class(SP500)
SPreturn <- SP500$r500
res <- jarque.bera.test(SPreturn)
names(res)
res$p.value

结果如下所示：

> class(SP500)
[1] "data.frame"
> library(Ecdat)
> data(SP500)
> class(SP500)
[1] "data.frame"
> SPreturn <- SP500$r500
> res <- jarque.bera.test(SPreturn)
> names(res)
[1] "statistic" "parameter" "p.value"   "method"    "data.name"
> res$p.value
[1] 0

正态分布检验补充案例

words1 = rnorm(100)
words2 = rt(100, df=3)

plot(density(words1))
lines(density(words2))
shapiro.test(words1)
shapiro.test(words2)

par(mfrow=c(2,1))
qqnorm(words1);qqline(words1, col = 2)
qqnorm(words2);qqline(words2, col = 2)

words1与words2密度图如下所示：

从图上大致可以看出来words1与标准正态图比较接近。

用shapiro检验后的结果如下所示：

从计算结果来看，words1是p值为0.09861，大于0.05，不能否定原假定，因此words1符合正态分布。同理,words2不符合正态分布。

QQ图的结果如下所示：

从QQ图中可以看出，words1符合正态分布，而words2不符合正态分布。

参考资料

1. 每日一生信—富集分析（超几何分布）（Fisher’s Exact Test）
2. R 超几何分布
3. Wikipedia:Hypergeometric distribution
4. DawnGriffiths. 深入浅出统计学[M]. 电子工业出版社, 2012.
5. 高桥信陈刚. 漫画统计学[M]. 科学出版社, 2009.
6. 拉里·戈尼克, 沃尔科特·史密斯. 漫画统计学入门[M]. 辽宁教育出版社, 2002.
7. 陈希孺.概率论与数理统计[M].中国科学技术大学出版社,1992.