如何得到服从正态分布的样本

之前介绍了计算机如何产生随机数，同时也整理了常用的概率分布。现在我们来考虑这样的问题，如何得到服从正态分布的样本？，在计算机如何产生随机数中我们用线性同余发生器生成了随机数，这些随机数实际上就是服从 $U(0,1)$ 的均匀分布，意味在 $(0,1)$ 范围内，每个样本被抽取概率是一样的。

对于正态分布的样本来说，每个样本被抽取的概率是不同的。根据正态分布的曲线我们知道，在平均值附近的点被抽取的概率是最大的，越往两边走，被抽取的样本的概率越低，于是形成了钟形曲线。

采样

采样的方式具有较强的通用性，一般用于比较复杂的概率分布。常用的采样方法有接受-拒绝采样，大名鼎鼎的马尔可夫蒙特卡洛采样（MCMC采样），其典型的代表是Metropolis-Hasting采样，Gibbs采样。这里我们先挖个坑留着，后续继续介绍。因为对于正态分布的采样虽然也能用这些采样的方法，但是效率不是最高的。用复杂的模型处理简单的问题，不太合理哈。

中心极限定理生成

中心极限定理描述的是多个独立的随机变量求和近似服从正态分布。之前我们得到了均匀分布，可以尝试独立地从均匀分布中采样多个变量，求和以后可以得到近似服从正态分布的样本。中心极限定理有多种描述方式，下面仅介绍其中一种。

独立同分布中心极限定理， $X_1,X_2,\cdots ,X_n$ ,为独立同分布的随机变量序列，均值 $\mu$ 为方差 $\sigma^2$ 均存在。

Z_n = \frac{\sum_i^n X_i - n\mu}{\sigma \sqrt n}

\lim_{n\rightarrow \infty} P(Z_n\le x) = \int_{-\infty} ^ x \frac{1}{\sqrt{2\pi}} \mathrm e^{-\frac{x^2}{2}} \mathrm dx

意味着 $Z_n$ 近似服从标准正态分布，随机变量的个数越多，越近似服从正态分布。高尔顿板就是典型的中心极限定理的应用，每个小球独立的服从二项分布，小球每次碰到杆随机的选择一次方向，最终的结果近似服从正态分布。

接下来用程序来生成正态分布的样本

import numpy as np
from matplotlib import pyplot as plt

mu = 0.5  # 期望
s = 1 / 12  # 方差
sample_size = 10000  # 样本数
n = 1000  # 随机变量个数

total = np.sum(np.random.random((n, sample_size)), axis=0)  # 随机变量求和
normal_dist = (total - n * mu) / np.sqrt(s * n)  # 近似服从正态分布
plt.figure(0)
plt.hist(normal_dist, bins=np.linspace(-4, 4, 100), facecolor='green')
plt.show()

可以看到 $n=1$ 时其实就是均匀分布，随着 $n$ 逐渐增大，直方图轮廓越来越接近正态分布了。但是这样生成正态分布样本的效率是很低的，因为我们需要生成大量的随机变量，每个随机变量也要采样大量的样本，才能得到正态分布的样本。

逆变换法生成

我们已经有了均匀分布的样本，考虑建立映射把均匀分布的样本投影到正态分布。这样就实现了从正态分布采样了，这种方法称为Box–Muller变换，在实际中应用很多，例如Python的源码random模块生成正态分布的随机数就是用的这个方法。

逆变换定理：设 $X$ 为连续型随机变量，取值于区间 $(a,b)$ (可包括 $\infty$ 和端点), $X$ 的概率密度在 $(a,b)$ 上取正值， $F(X)$ 的分布函数为 $F(X),U \sim \text{U}(0,1)$ 则 $Y = F^{-1}(U) \sim F(\cdot)$

假定 $X、Y$ 服从均值为0，方差为1的高斯分布，且相互独立，它们的概率密度如下。

p(X) = \frac 1 {\sqrt {2 \pi}} \mathrm e^{-\frac {X^2} 2}

p(Y) = \frac 1 {\sqrt {2 \pi}} \mathrm e^{-\frac {Y^2} 2}

P(X,Y) = P(X)P(Y) = \frac 1 {2 \pi} \mathrm e^{-\frac {X^2+Y^2} 2}

然后根据极坐标变换 $X= R\cos\theta, Y = R\sin\theta$ 则概率密度变成:

p(R, \theta) = \frac 1 {2\pi} \mathrm e^{-\frac {R^2} 2}R

然后分别得到 $r, \theta$ 的累积分布函数

F_\phi (\phi) = P( \theta <= \phi) = \int_0^{\phi} \int_0^{+\infty} \frac 1 {2\pi} e^{-\frac {R^2} 2}R \ \mathrm dR \mathrm d\theta = \frac \phi {2 \pi}

F_r( R \leq r) = \int_0^r \int_0^{2\pi} \frac 1 {2\pi} e^{-\frac {R^2} 2}R\ \mathrm d\theta \mathrm dR =1- \mathrm e^{-\frac {r^2} 2}

根据逆变换定理

R = F^{-1}(U_1) = \sqrt{-2\ln(1-U_1)}=\sqrt{-2\ln U_1}

\phi = 2\pi U_2

上面 $U_1, U_2\sim U(0,1)$ ,带入到 $X,Y$ 即可得到服从正态分布的样本

X = \sqrt{-2\ln U_1} \cos (2\pi U_2)

Y = \sqrt{-2\ln U_1} \sin (2\pi U_2)

接下来验证一下生成样本的效果。

import numpy as np
from matplotlib import pyplot as plt


def get_normal_dist(sample_size):
    u1 = np.random.uniform(size=sample_size)
    u2 = np.random.uniform(size=sample_size)

    return np.sqrt(-2 * np.log(u1)) * np.cos(2 * np.pi * u2)


plt.hist(get_normal_dist(10000),
         bins=np.linspace(-4, 4, 100),
         facecolor='green')
plt.show()

总结

本文主要介绍生成正态分布样本的方法，通用的采样方法，利用中心极限定理生成，逆变换法生成，其中重点介绍了后两种方法，并且给出了程序模拟。我做公众号的目的也是为了介绍数学编程，其实编程很大程度上就是数学的应用，把书本上的公式用在实际中，是不是能感受到数学与编程的魅力。

参考

https://cosx.org/2015/06/generating-normal-distr-variates

土豆不吃鱼