采样（二）：从正态分布采样

计算机程序的运行是确定性的，即每一步都有一个明确的描述，如何在确定性下生成随机的内容？这似乎是个自相矛盾的问题。

事实上，这个怀疑是对的，在没有外部不确定性的输入下，计算机程序永远以确定性的方式运行。为让计算机程序稳定地生成随机内容，需要稳定的不确定性的输入。而这样做的成本非常高。于是，为何不抛开这个直观思路？让计算机生成一些看起来随机的确定性序列。易知，伪随机数有一个非常好的优点，可重复。比如我们再训练神经网络时，要初始化大量参数，伪随机数则可以让我们获得同样的初始化参数，从而避免对此训练。

采样

采样，可以看做是从给定分布中采集样本，也可以看做是从给定数据集，通过某种策略，采集构造新的数据集的方法。本文讲述从均匀分布采样、从指数分布采样、从正态分布、Laplace分布采样等常见分布采样方法。

首先我们谈谈分布变换，即有已知分布的样本集，通过一定的方法使其变换为符合其他分布的样本集。比如，我现在有一个均匀分布生成器，获得一系列样本$x_1, x_2, \dots, x_n \sim U[0, 1]$，那么通过$g(x)$，它如何变换为其他的分布呢？比如变换为正态分布，

$y_i = g(x_i) \sim N(0, 1)$

这就是分布变换的问题。一般来说，简单的分布容易采样，在获得简单分布的样本后，通过分布变换获得复杂分布的样本。但是对于一般的分布函数，采样方法都不是普遍适用的，因此需要具体分布具体讨论，为每种分布设计最高效的采样方法。

逆变换采样

累积分布函数（CDF）是个单调函数，那么累积分布函数的反函数为，

$F^{-1}(u) = \inf \left \{ x |F(u) \geq x \right \}$

如果$u \sim U[0, 1]$，那么$F^{-1}(u)$服从累积分布$F(x)$。

证明，

$\begin{align} P(F^{-1}(u) \le x) &= P(u \le F(x)) \newline &= P(u \le y) \Big|_{y = F(x)} \newline &= y|_{y = F(x)} \newline &= F(x) \end{align}$

这里有一个细节需要说明，由于$u$采样自均匀分布，因此容易知道$P(u \le y) = y$。这一点是以上推导的关键细节。对于很多分布来说，逆函数$F^{-1}(u)$并不容易计算，因此很多情况下无法直接使用逆变换采样。

注意到，

$F(F^{-1}(u)) = u$

也就是说，如果获得累积分布为$F(x)$的随机变量$t$，那么$F(t)$为均匀分布。例如$x_i \sim N(0, 1)$，那么标准正态分布的累积分布函数（CDF）$\Phi(x_i)$服从均匀分布。不过$\Phi(x)$无法直接计算，可以使用其近似形式，

$\Phi(x) \approx \frac{1}{2} \Big( 1 + \tanh \left[ \sqrt{\frac{2}{\pi}} (x + 0.044715x^3)\right] \Big)$

于是，对于$x \sim N(0, 1)$代入上式即可获得$U[0, 1]$样本。

正态分布

正态分布的概率密度函数，

$f(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}e^{-(x-\mu)^{2}/2\sigma^{2}}$

均值和方差分别为$\mu$和$\sigma^{2}$。

取 $z=\frac{x-\mu}{\sigma}$ 化成标准正态分布，

$f(x)=\frac{1}{\sqrt{2 \pi}} \mathrm{e}^{-x^{2} / 2}$

如果从标准正态分布采样$u_i$，那么$\mu + u_i \times \sigma$就是均值和方差分别为$\mu$和$\sigma^{2}$的随机变量。

逆变换采样法

计算累计分布函数，

$\begin{align} \Phi(x) &= \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}}e^{-x^2/2}dx \newline &= \int_{-\infty}^{0} \frac{1}{\sqrt{2\pi}}e^{-x^2/2}dx + \int_{0}^{x} \frac{1}{\sqrt{2\pi}}e^{-x^2/2}dx \newline &= \frac{1}{2} + \frac{1}{2}\operatorname{erf}(\frac{x}{\sqrt{2}}) \end{align}$

对于一般正态分布也容易得到

$F(x)=\Phi \left({\frac {x-\mu }{\sigma }}\right)={\frac {1}{2}}\left[1+\operatorname {erf} \left({\frac {x-\mu }{\sigma {\sqrt {2}}}}\right)\right]$

然后你会发现其逆$F^{-1}(x)$没那么容易求。scipy这个科学计算库给我们提供了一种实现。Python实践如下，

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# 正态分布采样的逆变换法，函数为scipy.stats.norm.pdf
# 原理为分段近似

def gen_normal_curve(rs, u, s):
    x = np.linspace(np.min(rs), np.max(rs), len(rs))
    rs = 1 / (np.sqrt(2*np.pi)*s) * np.exp(-(x-u)**2/(2*s**2))
    return x, rs

x = stats.norm.ppf(np.random.uniform(size=10000))
print(stats.normaltest(x))
plt.hist(x, bins=100, density=True)
mu = 0
sigma = 1
x, c = gen_normal_curve(x, u=mu, s=sigma)
plt.plot(x, c, color="red", label="$N({},{}^2)$".format(mu, sigma))
plt.legend(loc="upper left")
plt.show()

可视化如下，

可以看到采样样本（蓝色部分）和标准正态分布拟合得非常好。

另外一种思路是寻找一个可逆的函数逼近误差函数。注意到$\operatorname{erf}(x)$可以用可逆的函数光滑逼近，

$\begin{align} \Phi(x) &= \frac{1}{2} + \frac{1}{2}\operatorname{erf}(\frac{x}{\sqrt{2}}) \newline &\approx \frac{1}{1 + e^{-2kx}} \end{align}$

求逆，解得$x$为，

$x = -\frac{1}{2k} \ln( \frac{1}{u} - 1)$

可视化这种采样效果，

可以看到，还是有较大差异。

此外$\operatorname{erf}(x)$还有一个不错的光滑逼近，

$\operatorname{erf}(\frac{x}{\sqrt{2}}) \approx \tanh \left[ \sqrt{\frac{2}{\pi}} (x + 0.044715x^3)\right]$

这个出发点，然后求$\tanh$的逆，然后解一元三次方程，获得逼近函数的逆，作为$\operatorname{erf}(x)$的逆的近似。

正态分布的CDF无法直接求解反变换函数，但是可以通过一定的技巧绕过这个问题。例如使用 ziggurat 有效地求累积分布函数的逆函数。

Box-Muller 方法

既然$\operatorname{erf}(x)$和$\Phi(x)$的逆变换难求，那么我们对其升维后再试试使用逆变换的思路。这是机器学习的思路：低维无法解决问题，升维来解决。

这里的具体思路是$\Phi(x)$难求逆函数，但是二维$\Phi(x)\Phi(y)$容易求逆函数，获得逆函数后通过逆变换采样得到正态分布的采样公式，这种方法称为Box-Muller方法。

二维的标准正态分布为，

$p(x,y)=p(x)p(y) = \frac{1}{2 \pi} \mathrm{e}^{-\left(x^{2}+y^{2}\right) / 2}$

其累积分布函数通过极坐标变换容易求解。这里另，

$x = r\cos(\theta) \newline y = r \sin(\theta)$

于是，我们解决$r$与$\cos(\theta),\sin(\theta)$的采样问题即可解决正态分布采样问题。

考虑如下积分计算，

$\begin{align} \Phi(x) \Phi(y) &= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} e^{\frac{-x^2}{2}} dx \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{\frac{-y^2}{2}} dy \newline &= \frac{1}{2\pi} \int_{r=0}^{\infty} \int_{\theta=0}^{2\pi} e^{\frac{-r^2}{2}} r dr d\theta \newline &=\frac{1}{2\pi} \int_{\theta=0}^{2\pi} d\theta \int_{r=0}^{\infty} e^{\frac{-r^2}{2}} r dr \newline &= \int_{r=0}^{\infty} e^{\frac{-r^2}{2}} r dr\newline &= 1 \end{align}$

该积分值为1，整体上可以看做是某个随机变量的概率密度函数。根据这个结果，我们使用逆变换思路。对于随机变量$R$，根据以上结果有，

$P(R \le r) = \int_{x=0}^{r} e^{\frac{-x^2}{2}} x dx\ = 1 - e^{-\frac{r^2}{2}} = u_1 \sim U[0, 1]$

对于指数分布有，

$P(x)=\operatorname{P}(X \leqslant x)=\int_{0}^{x} p(x) \mathrm{d} x=1-\mathrm{e}^{-a x}$

于是，这里把 $r^{2}$ 看作一个整体，$P(r)$ 实质是参数$a = \frac{1}{2}$的指数分布，于是逆变换采样思路有，

$r^2 = -2 \ln(1 - u_1)$

也就有$r$的采样方法，

$r = \sqrt{-2 \ln(1 - u_1)}$

考虑到$\theta$的取值范围为$2\pi u_2$，于是根据$x = r\cos(\theta), y = r \sin(\theta)$，有正态分布采样方法，

$\begin{array} {l}{n_{1}=r \cos 2 \pi u_{2}=\sqrt{-2 \ln \left(1-u_{1}\right)} \cos 2 \pi u_{2}} \newline {n_{2}=r \sin 2 \pi u_{2}=\sqrt{-2 \ln \left(1-u_{1}\right)} \sin 2 \pi u_{2}} \end{array}$

生成正态分布任取其中一条公式即可。Box-Muller 方法需要进行三角运算，计算量较大。通过一定的技巧可以把其去掉三角运算。

以上技巧如何逆着来，即可通过正态分布生成器获得均匀分布生成器，$n_1,n_2$独立采样自标准正态分布，那么

$u = 1 - e^{-\frac{n_1^2 + n_2 ^2}{2}}$

服从均匀分布$U[0, 1]$。

修正的 Box-Muller 方法

假设$x \sim U[-1, 1], y \sim U[-1, 1]$，即在矩形上采样$(x, y)$，对于满足$r = x^{2} + y^{2} \le 1$的样本留下，否则丢弃，这样获得的样本$r$满足均匀分布，因为这相当于等可能地落到圆内。这个方法也称为Marsaglia polar method。

于是，$\frac{x}{\sqrt{r}}$可以作为$\cos 2 \pi u_{2}$的替代；$\frac{y}{\sqrt{r}}$可以作为$\sin 2 \pi u_{2}$的替代。

因此有，

$n_{1}= \sqrt{-2 \ln \left(1-u_{1}\right)} \cos 2 \pi u_{2} = \frac{x}{\sqrt{r}} \times \sqrt{-2\ln r} \newline n_{2}= \sqrt{-2 \ln \left(1-u_{1}\right)} \sin 2 \pi u_{2} = \frac{y}{\sqrt{r}} \times \sqrt{-2\ln r}$

修正的 Box-Muller 方法有一个好处是避免进行三角运算，提升计算效率。

中心极限定理思路

中心极限定理表明多个独立统计量（可以具有不同分布）的平均值满足正态分布。以均匀分布为例，这里均匀分布$x_{i} \sim U[0, 1]$，易知均值为$\frac{1}{2}$，方差为$\frac{1}{12}$，从中采样一系列样本$x_1, \dots, x_n$，有

$\begin{align} y &= \frac{\displaystyle \sum_{i=1}^{n} x_i - E\Big[\displaystyle \sum_{i=1}^{n} x_i \Big]}{\sqrt{\operatorname{Var}(\displaystyle \sum_{i=1}^{n} x_i) }} \newline &= \frac{\displaystyle \sum_{i=1}^{n} x_i - \frac{n}{2}}{\sqrt{\frac{n}{12}}} \newline &\sim N(0, 1^2) \end{align}$

也就有，

$\frac{1}{n} \sum_{i=1}^{n} x_{i} \sim N(\frac{n}{2}, \frac{n}{12})$

以上是对于采样分布均为均匀分布的特殊情况。中心极限定理表明多个具有不同分布独立统计量也满足条件，这里就不进行数学推导了，Python实践编程验证一下。

一下的实现从正态分布、泊松分布、拉普拉斯分布、指数分布等中采样样本，

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

poisson = lambda: np.random.poisson(lam=100)
uniform = lambda: np.random.uniform()
normal = lambda: np.random.normal()
laplace = lambda: np.random.laplace()
gumbel = lambda: np.random.gumbel()
exponential = lambda: np.random.exponential()

probs = [poisson, uniform, normal, laplace, gumbel, exponential]

def gen_random(n=100, p=None):
    nums = []
    for _ in range(n):
        f = np.random.choice(probs, p=p)
        nums.append(f())
    return np.mean(nums)

size = 10000
nums = [gen_random(100) for _ in range(size)]
plt.hist(nums, bins=100, density=True)

mu = np.mean(nums)
sigma = np.std(nums)
x = np.linspace(mu - 4*sigma, mu + 4*sigma, size)
y = stats.norm.pdf(x, mu, sigma)
plt.plot(x, y, 'r')
plt.show()