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

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

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

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

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

yi=g(xi)∼N(0,1)yi=g(xi)∼N(0,1)

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

逆变换采样

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

F−1(u)=inf{x|F(u)≥x}F−1(u)=inf{x|F(u)≥x}

如果u∼U[0,1]u∼U[0,1]，那么F−1(u)F−1(u)服从累积分布F(x)F(x)。

P(F−1(u)≤x)=P(u≤F(x))=P(u≤y)∣∣y=F(x)=y|y=F(x)=F(x)P(F−1(u)≤x)=P(u≤F(x))=P(u≤y)|y=F(x)=y|y=F(x)=F(x)

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

F(F−1(u))=uF(F−1(u))=u

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

Φ(x)≈12(1+tanh[√2π(x+0.044715x3)])Φ(x)≈12(1+tanh⁡[2π(x+0.044715x3)])

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

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

f(x|μ,σ2)=1√2πσe−(x−μ)2/2σ2f(x|μ,σ2)=12πσe−(x−μ)2/2σ2

均值和方差分别为μμ和σ2σ2。

取 z=x−μσz=x−μσ​ 化成标准正态分布，

f(x)=1√2πe−x2/2f(x)=12πe−x2/2

如果从标准正态分布采样uiui，那么μ+ui×σμ+ui×σ就是均值和方差分别为μμ和σ2σ2的随机变量。

逆变换采样法

计算累计分布函数，

Φ(x)=∫x−∞1√2πe−x2/2dx=∫0−∞1√2πe−x2/2dx+∫x01√2πe−x2/2dx=12+12erf(x√2)Φ(x)=∫−∞x12πe−x2/2dx=∫−∞012πe−x2/2dx+∫0x12πe−x2/2dx=12+12erf⁡(x2)

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

F(x)=Φ(x−μσ)=12[1+erf(x−μσ√2)]F(x)=Φ(x−μσ)=12[1+erf⁡(x−μσ2)]

然后你会发现其逆F−1(x)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()

可视化如下，

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

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

Φ(x)=12+12erf(x√2)≈11+e−2kxΦ(x)=12+12erf⁡(x2)≈11+e−2kx

求逆，解得xx​为，

x=−12kln(1u−1)x=−12kln⁡(1u−1)

可视化这种采样效果，

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

此外erf(x)erf⁡(x)还有一个不错的光滑逼近，

erf(x√2)≈tanh[√2π(x+0.044715x3)]erf⁡(x2)≈tanh⁡[2π(x+0.044715x3)]

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

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

Box-Muller 方法

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

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

二维的标准正态分布为，

p(x,y)=p(x)p(y)=12πe−(x2+y2)/2p(x,y)=p(x)p(y)=12πe−(x2+y2)/2

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

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

于是，我们解决rr与cos(θ),sin(θ)cos⁡(θ),sin⁡(θ)的采样问题即可解决正态分布采样问题。

考虑如下积分计算，

Φ(x)Φ(y)=1√2π∫∞−∞e−x22dx1√2π∫∞−∞e−y22dy=12π∫∞r=0∫2πθ=0e−r22rdrdθ=12π∫2πθ=0dθ∫∞r=0e−r22rdr=∫∞r=0e−r22rdr=1Φ(x)Φ(y)=12π∫−∞∞e−x22dx12π∫−∞∞e−y22dy=12π∫r=0∞∫θ=02πe−r22rdrdθ=12π∫θ=02πdθ∫r=0∞e−r22rdr=∫r=0∞e−r22rdr=1

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

P(R≤r)=∫rx=0e−x22xdx =1−e−r22=u1∼U[0,1]P(R≤r)=∫x=0re−x22xdx=1−e−r22=u1∼U[0,1]

对于指数分布有，

P(x)=P(X⩽x)=∫x0p(x)dx=1−e−axP(x)=P⁡(X⩽x)=∫0xp(x)dx=1−e−ax

于是，这里把 r2r2 看作一个整体，P(r)P(r) 实质是参数a=12a=12的指数分布，于是逆变换采样思路有，

r2=−2ln(1−u1)r2=−2ln⁡(1−u1)

也就有rr的采样方法，

r=√−2ln(1−u1)r=−2ln⁡(1−u1)

考虑到θθ​的取值范围为2πu22πu2​，于是根据x=rcos(θ),y=rsin(θ)x=rcos⁡(θ),y=rsin⁡(θ)​，有正态分布采样方法，

n1=rcos2πu2=√−2ln(1−u1)cos2πu2n2=rsin2πu2=√−2ln(1−u1)sin2πu2n1=rcos⁡2πu2=−2ln⁡(1−u1)cos⁡2πu2n2=rsin⁡2πu2=−2ln⁡(1−u1)sin⁡2πu2

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

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

u=1−e−n21+n222u=1−e−n12+n222

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

修正的 Box-Muller 方法

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

于是，x√rxr可以作为cos2πu2cos⁡2πu2的替代；y√ryr可以作为sin2πu2sin⁡2πu2​的替代。

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}n_{1}= \sqrt{-2 \ln \left(1-u_{1}\right)} \cos 2 \pi u_{2} = \frac{x}{\sqrt{r}} \times \sqrt{-2\ln r} \newlinen_{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 方法有一个好处是避免进行三角运算，提升计算效率。

中心极限定理思路

中心极限定理表明多个独立统计量（可以具有不同分布）的平均值满足正态分布。以均匀分布为例，这里均匀分布xi∼U[0,1]xi∼U[0,1]​，易知均值为1212​，方差为112112​，从中采样一系列样本x1,…,xnx1,…,xn，有

y=n∑i=1xi−E[n∑i=1xi]√Var(n∑i=1xi)=n∑i=1xi−n2√n12∼N(0,12)y=∑i=1nxi−E[∑i=1nxi]Var⁡(∑i=1nxi)=∑i=1nxi−n2n12∼N(0,12) 1nn∑i=1xi∼N(n2,n12)1n∑i=1nxi∼N(n2,n12)

以上是对于采样分布均为均匀分布的特殊情况。中心极限定理表明多个具有不同分布独立统计量也满足条件，这里就不进行数学推导了，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()

可视化如下，

可以看到采样样本（蓝色部分）和正态分布（红色曲线）拟合得较好。

这种方法不建议工程环境使用：

• 这个思路更多是理论意义，实际计算效率不高
• 正态检验pvalue非常小，不显著

实现更新到：sampling-from-distribution，可能会根据需要持续更新~

本文介绍了一些正态分布采样方法，具体如下：

• Box-Muller
• 修正的Box-Muller
• 接受-拒绝采样
• 累积分布函数的逆或近似逆

此外，还有accept-reject抽样，可以用于复杂形式的概率密度函数的抽样。

[1] https://en.wikipedia.org/wiki/Lehmer_random_number_generator

[2] https://en.wikipedia.org/wiki/Inverse_transform_sampling

[3] https://en.wikipedia.org/wiki/Diehard_tests

转载请包括本文地址：https://allenwind.github.io/blog/10395
更多文章请参考：https://allenwind.github.io/blog/archives/