在学习了一些基本的统计变量生成法之后，这次我们来看看如何生成正态分布。它就是大名鼎鼎的 Box-Muller 方法，Box-Muller 的理解过程可以体会到统计模拟的一些精妙思想。

• 从零构建统计随机变量生成器之 离散基础篇
• 从零构建统计随机变量生成器之 用逆变换采样方法构建随机变量生成器
• 深入 LeetCode 470 拒绝采样，状态转移图求期望和一道经典统计求期望题目
• 从零构建统计随机变量生成器之 正态分布 Box-Muller方法

尝试逆变换方法

我们先尝试通过标准的逆变换方法来生成正态分布。

正态分布的 PDF 表达式为

对应的函数图形是钟形曲线

根据 PDF，其 CDF 的积分形式为 \[ \Phi(x)=\frac{1}{\sqrt{2 \pi}} \int_{-\infty}^{x} e^{-t^{2} / 2} d t \]

和所有 PDF CDF 关系一样，\(\Phi(x)\) 表示 \(f_Z\) 累积到 \(x\) 点的面积。

很不幸的是，\(\Phi(x)\) 无法写出一般数学表达式，因而也无法直接用逆变换方法。

二维映射到一维

我们知道，高维正态分布有特殊的性质：它的每一维的分量都是正态分布；单个维度对于其他维度的条件概率分布也是正态分布。

用图来理解这两条性质就是，对于下图的二维正态分布 $ x = [x_1, x_2]^T $，单独的 \(x_1\) 和 \(x_2\) 都服从一维正态分布。

条件概率 \(p(x_2|x_1 \approx1)\) 的PDF 对应图中的红线，显然也是一维正态分布。

写一段简单的代码验证二维正态分布的单个分量服从正态分布。

代码中，我们用np.random.normal生成了 10000 个服从二维正态分布的 x, y 点，然后我们丢弃 y，只保留 x，并画出 10000 个 x 的分布。

1
2
3
4
5
6
7
8

def plot_normal_1d():
    x, _ = np.random.normal(loc=0, scale=1, size=(2, 10000))
    import seaborn as sns
    sns.distplot(x, hist=True, kde=True, bins=100, color='darkblue',
                 hist_kws={'edgecolor': 'black'},
                 kde_kws={'linewidth': 4})
    plt.title('PDF Normal 1D from 2D')
    plt.show()

Box-Muller 原理

虽然无法直接用逆变换方法生成一维正态分布，但我们却能通过先生成二维的正态分布，利用上面一节的性质，生成一维正态分布。

而 Box-Muller 就是巧妙生成二维正态分布样本点的方法。

首先，我们来看看二维正态分布可以认为是两个维度是独立的，每个维度都是正态分布。此时，其 PDF 可以写成两个一维正态分布 PDF 的乘积。

这种写法表明，二维正态分布仅用一个 r 向量就可以充分表达。注意，r 是向量，不仅有大小还有角度，有两个分量。这两个分量本质上是独立的，这就是 Box-Muller 方法的巧妙之处。也就是，Box-Muller 通过角度和半径大小两个分量的独立性分别单独生成并转换成 (x, y) 对。

角度分量是在 \(2\pi\) 范围均匀采样，这一点比较直觉好理解。

再来看看半径分量 r。我们令 \[ s = {r^2 \over 2} \Longrightarrow r = \sqrt{2s} \]

则 s 服从指数分布 \(\lambda=1\) 。

不信么？我们不妨来做个模拟实验，下图是模拟 10000次二维正态分布 (x, y) 点后转换成 s 的分布。

模拟和plot 代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

def plot_r_squared():
    def gen_normal_samples(n):
        x, y = np.random.normal(loc=0, scale=1, size=(2, n))
        return x, y

    x, y = gen_normal_samples(10000)
    s = (x * x + y * y)/2
    plot_dist_1d(s, title='PDF $s = {{x^2 + y^2}\over{2}} \sim exp(1)$')
    

def plot_dist_1d(X, title='PDF '):
    import seaborn as sns
    plt.rcParams.update({
        "text.usetex": True,
        "font.family": "sans-serif",
        "font.sans-serif": ["Helvetica"]})
    sns.distplot(X, hist=True, kde=True, bins=100, color='darkblue',
                 hist_kws={'edgecolor': 'black'},
                 kde_kws={'linewidth': 4})
    plt.title(title)
    plt.show()    

确信了 s 符合指数分布，根据指数分布的 PDF，可以推出二维正态 PDF中的 $ e^{-r2/2}$ 也符合指数分布，即 \[ s \sim \exp(1) \Longrightarrow e^{-r^2/2} \sim \exp(1) \]

至此，总结一下Box-Muller方法。我们视二维正态分布PDF为独立两部分的乘积，第一部分是在 \(2\pi\) 范围中的均匀分布，代表了二维平面中的角度 \(\theta\)，第二部分为 \(\lambda=1\) 的指数分布，代表半径大小。

Box-Muller 方法通过两个服从 [0, 1] 均匀分布的样本 u1和u2，转换成独立的角度和半径样本，具体过程如下

1. 生成 [0, 1] 的均匀分布 u1，利用逆变换采样方法转换成 exp(1) 样本，此为二维平面点半径 r

2. 生成 [0, 1] 的均匀分布 u2，乘以 \(2\pi\)，即为样本点的角度 \(\theta\)

3. 将 r 和 \(\theta\) 转换成 x, y 坐标下的点。

理解了整个过程的意义，下面的代码就很直白。

1
2
3
4
5
6
7
8
9
10

def normal_box_muller():
    import random
    from math import sqrt, log, pi, cos, sin
    u1 = random.random()
    u2 = random.random()
    r = sqrt(-2 * log(u1))
    theta = 2 * pi * u2
    z0 = r * cos(theta)
    z1 = r * sin(theta)
    return z0, z1

接下来，我们来看看 Box-Muller 法生成的二维标准正态分布动画吧

拒绝采样极坐标方法

Box-Muller 方法还有一种形式，称为极坐标形式，属于拒绝采样方法。

1. 生成独立的 u, v 和 s

分别生成 [0, 1] 均匀分布 u 和 v。令 \(s = r^2 = u^2 + v^2\)。如果 s = 0或 s ≥ 1，则丢弃 u 和 v ，并尝试另一对 (u , v)。因为 u 和 v 是均匀分布的，并且因为只允许单位圆内的点，所以 s 的值也将均匀分布在开区间 (0, 1) 中。注意，这里的 s 的意义虽然也为半径，但不同于基本方法中的 s。这里 s 取值范围为 (0, 1) ，目的是通过 s 生成指数分布，而基本方法中的 s 取值范围为 [0, +∞]，表示二维正态分布 PDF 采样点的半径。复用符号 s 的原因是为了对应维基百科中关于基本方法和极坐标方法的数学描述。

我们用代码来验证 s 服从 (0, 1) 范围上的均匀分布。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

def gen_polar_s():
    import random
    while True:
        u = random.uniform(-1, 1)
        v = random.uniform(-1, 1)
        s = u * u + v * v
        if s >= 1.0 or s == 0.0:
            continue
        return s


def plot_polar_s():
    s = [gen_polar_s() for _ in range(10000) ]
    plot_dist_1d(s, title='PDF Polar $s = u^2 + v^2$')

2. 将 u, v, s 转换成 x, y

若将 $s = R^2 uniform(0, 1) $ 看成是基本方法中的 u1，就可以用同样的方式转换成指数分布，用以代表二维PDF的半径。

同时，根据下图，\(\cos \theta\) 和 \(\sin \theta\) 可以直接用 u, v, R 表示出来，并不需要通过三角函数显示计算出 \(\theta\)。有了半径， \(\cos \theta\) 和 \(\sin \theta\) ，则可以直接计算出 x, y 坐标，（下面用 \(z_0, z_1\) 代替 \(x, y\)）。

\[ z_{0}=\sqrt{-2 \ln U_{1}} \cos \left(2 \pi U_{2}\right)=\sqrt{-2 \ln s}\left(\frac{u}{\sqrt{s}}\right)=u \cdot \sqrt{\frac{-2 \ln s}{s}} \]

\[ z_{1}=\sqrt{-2 \ln U_{1}} \sin \left(2 \pi U_{2}\right)=\sqrt{-2 \ln s}\left(\frac{v}{\sqrt{s}}\right)=v \cdot \sqrt{\frac{-2 \ln s}{s}} \]

同样，Box-Muller 极坐标方法的代码和公式一致。

1
2
3
4
5
6
7
8
9
10
11
12

def normal_box_muller_polar():
    import random
    from math import sqrt, log
    while True:
        u = random.uniform(-1, 1)
        v = random.uniform(-1, 1)
        s = u * u + v * v
        if s >= 1.0 or s == 0.0:
            continue
        z0 = u * sqrt(-2 * log(s) / s)
        z1 = v * sqrt(-2 * log(s) / s)
        return z0, z1

拒绝采样的效率

极坐标方法与基本方法的不同之处在于它是一种拒绝采样。因此，它会丢弃一些生成的随机数，但可能比基本方法更快，因为它计算更简单：避免使用昂贵的三角函数，并且在数值上更稳健。极坐标方法丢弃了生成总输入对的 1 − π /4 ≈ 21.46%，即需要 4/ π ≈ 1.2732 个输入随机数，输出一个随机采样。

导读：极大似然估计(MLE) 是统计机器学习中最基本的概念，但是能真正全面深入地理解它的性质和背后和其他基本理论的关系不是件容易的事情。极大似然估计和以下概念都有着紧密的联系：随机变量，无偏性质（unbiasedness），一致估计（consistent），asymptotic normality，最优化（optimization），Fisher Information，MAP（最大后验估计），KL-Divergence，sufficient statistics等。在众多阐述 MLE 的文章或者课程中，总体来说都比较抽象，注重公式推导。本系列文章受3blue1brown 可视化教学的启发，坚持从第一性原理出发，通过数学原理结合模拟和动画，深入浅出地让读者理解极大似然估计。

相关链接

抛硬币问题

我们来思考这个老套问题，考虑手上有一枚硬币，旋转（抛）硬币得到正反面的概率固定（令正面概率为\(\theta^{\star}\)）但未知，我们如何能通过实验推测出 \(\theta^{\star}\)

朴素的想法是，不断尝试抛硬币，随着次数 n 的增多，正面的比例会趋近于 \(\theta^{\star}\)

对应到数学形式上，令我们对于 \(\theta^{\star}\) 的估计为 \(\hat{\theta}_{n}\)，则希望 \[ \hat{\theta}_n = {n_{head} \over n} \to \theta^{\star} \text{ as n } \to \infty \]

模拟试验代码

假设我们尝试了n次，每次的结果为 \(x_i\)，\(x_i\)为1（正面） 或 0（反面）。比如试了三次的结果是 [1, 0, 1]，则 \(x_1=1, x_2=0, x_3=1\)。一般，我们将观察到的数据写成向量形式

\[X=[x_1, x_2, x_3]^T=[1, 0, 1]^{T}\]

我们知道硬币的正反结果符合伯努利分布，也就是 \[ \begin{align*} P_{ber}(x;\theta) = \left\lbrace \begin{array}{r@{}l} \theta &\text{ if x=1} \\ 1-\theta &\text{ if x=0} \end{array} \right. \end{align*} \]

因为 x 只有0，1两种取值，因此上式也可以写成等价如下的不含条件分支的形式 \[ P_{ber} = \theta^x \cdot (1-\theta)^x \]

假设 \(\theta^{\star} = 0.7\)，如果做 n=10 次试验，结果应该比较接近7个1，3个0。

下面我们来模拟一下 n=10，看看结果如何。

下面代码的实现上我们直接使用了pytorch 内置的 bernoulli 函数生成 n 个随机变量实例

1
2
3
4
5

def gen_coins(theta, n=1):
    import torch
    theta_vec = torch.tensor(n*[theta])
    random_values = torch.bernoulli(theta_vec)
    return random_values

让我们来做三次 n=10 的试验

1
2
3
4
5
6

for i in range(3):
    coins = gen_coins(theta=0.7, n=10)
    print(f'trial {i}')
    print(f'head #: {sum(coins)}')
    print(f'tail #: {sum(1-coins)}')
    print()

能发现 7个1，3个0 确实是比较可能的结果。

1
2
3
4
5
6
7
8
9
10
11

trial 0
head #: 7.0
tail #: 3.0

trial 1
head #: 9.0
tail #: 1.0

trial 2
head #: 7.0
tail #: 3.0

生成概率

直觉告诉我们，当 \(\theta^{\star}=0.7\) 时，根据 $P_{ber}(x;) $，7个1，3个0 出现的概率应该是最大，6个1，4个0 或者 8个1，2个0 这两种情况出现概率稍小，其他的情况概率更小。通过基本概率和伯努利公式，重复 n 次试验 1和0出现的概率可以由下面公式算出。（注：7个1，3个0不是单一事件，需要乘以组合数算出实际概率）

\[ P_{X} = \theta^{heads} \cdot (1-\theta)^{tails} \cdot {n \choose heads} \]

画出图看的很明显，1出现7次的概率确实最大。

回到我们的问题，我们先假定 \(\theta^{\star} = 0.7\) 的硬币做 n=10 次试验的结果就是 7个1，3个0，或者具体序列为 [1, 0, 0, 1, 0, 1, 1, 1, 1, 1]。那么我们希望按照某种方法推测的估计值 \(\hat\theta\) 也为 0.7。

若将这个方法也记做 \(\hat\theta\)，它是\(X\) 的函数，即 \(\hat\theta(X=[1, 0, 0, 1, 0, 1, 1, 1, 1, 1]^T)=0.7\)

我们如何构建这个方法呢？很显然，\(X\) 中 1 的个数就可以胜任，\(\hat\theta=\bar X\)。这个方式确实是正确的，后面的文章我们也会证明它是MLE在伯努利分布参数估计时的计算方法。

但是伯努利分布参数估计的问题中是最简单的情况，背后对应的更一般的问题是：假设我们知道某个过程或者实验生成了某种分布 P，但是不知道它的参数 \(\theta\)，如何能通过反复的试验来推断 \(\theta\)，同时，我们希望随着试验次数的增多，\(\hat\theta\) 能逼近 \(\theta\)。

由于过程是有随机性，试验结果 \(X\) 并不能确定一定是从 \(\hat\theta\) 生成的，因此我们需要对所有 \(\theta\) 打分。对于抛硬币试验来说，我们穷举所有在 [0, 1] 范围内的 \(\theta\)，定义它的打分函数 \(f(\theta)\)，并且希望我们定义的 \(f(\theta;X=[1, 0, 0, 1, 0, 1, 1, 1, 1, 1]^T)\) 在 \(\theta=0.7\) 时得分最高。推广到一般场景，有如下性质 \[ f(\theta^\star;X) >= f(\theta;X) \]

如此，我们将推测参数问题转换成了优化问题 \[ \hat\theta = \theta^{\star} = \operatorname{argmax}_{\theta} f(\theta; X) = 0.7 \]

朴素方法

一种朴素的想法是，由于 \(\theta^\star=0.7\)，因此我们每次的结果应该稍微偏向 1，如果出现了 1，就记0.7分，出现了0，记0.3分，那么我们可以用10个结果的总分来定义总得分，即最大化函数

\[ \begin{equation*} \begin{aligned} &\operatorname{argmax}_{\theta} f(\theta) \\ =& \operatorname{argmax}_{\theta} P(x_1) + P(x_2) + ... + P(x_n) \\ =& \operatorname{argmax}_{\theta} P(x_1|\theta) + P(x_2|\theta) + ... + P(x_n|\theta) \\ =& \operatorname{argmax}_{\theta} \sum P(x_i|\theta) \\ \end{aligned} \end{equation*} \]

很可惜，我们定义的 f 并不符合 \(\theta=0.7\) 时取到最大的原则。下面画出了 \(\theta\) 在 [0, 1] 范围内 f 值，X 固定为 [1, 0, 0, 1, 0, 1, 1, 1, 1, 1]。显然，极值在 0.5 左右。

这种对于观察到的变量实例在整个参数空间打分的方法是最大似然方法的雏形。我们将每次试验结果对于不同 \(\theta\) 的打分就是似然函数的概念。

伯努利似然函数（Likelihood)

伯努利单个结果的似然函数 \(l(\theta)\) 视为 \(\theta\) 的函数，x视为给定值，它等价于概率质量函数 PMF

\[ l(\theta|x) = \theta^x \cdot (1-\theta)^x \]

极大似然估计(MLE)

有了单个结果的似然函数，我们如何定义 \(f(\theta)\) 呢？我们定义的 \(f(\theta)\) 需要满足，在 \(\theta^\star=0.7\) ，\(n=10\) 的情况下，试验最有可能的结果是 7 个1，3个0，此时 f 需要在 \(\theta=0.7\) 时取到最大值。

极大似然估计(MLE) 为我们定义了合理的 \(f(\theta)\) ，和朴素的想法类似，但是这次用单个结果的似然函数连乘而非连加 \[ L(\theta|X) = l(\theta|x_1) \cdot l(\theta|x_2) \cdot ...l(\theta|x_n) = \prod l(\theta|x_i) \]

我们再来看一下当 $X=[1, 0, 0, 1, 0, 1, 1, 1, 1, 1] $ 时 \(L\) 在 \(\theta\) 空间的取值情况，果然，MLE 能在 0.7时取到最大值。

对数似然函数

最大似然函数 $_{} L() $ 能让我们找到最可能的 \(\theta\)，但现实中，我们一般采用最大其 log 的形式。

\[ \begin{equation*} \begin{aligned} &\operatorname{argmax}_{\theta} \log L(\theta|X) \\ =& \operatorname{argmax}_{\theta} \log [l(\theta|x_1) \cdot l(\theta|x_2) \cdot ... \cdot l(\theta|x_n)] \\ =& \operatorname{argmax}_{\theta} \log l(\theta|x_1) + \log l(\theta|x_2) \cdot ... + \log l(\theta|x_n) \end{aligned} \end{equation*} \]

理论能证明，最大对数似然函数得到的极值等价于最大似然函数。但这么做有什么额外好处呢？

我们先将对数似然函数画出来

它的极大值也在 0.7，但是我们发现对数似然函数是个 concave 函数。在优化领域，最大化 concave 函数或者最小化 convex 函数可以有非常高效的解法。再仔细看之前的似然函数，它并不是一个 concave 函数。另一个非常重要的好处是，随着 n 的增大，连乘会导致浮点数 underflow，而单个点的对数似然函数的和的形式就不会有这个问题。

Pytorch MLE 代码

就让我们来实践一下，通过 pytorch 梯度上升来找到极值点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

from stats.coin import gen_coins
from collections import deque


def train(num_head: int, num_tail: int) -> float:
    import torch
    theta = torch.tensor(0.5, requires_grad=True)

    recent = deque(3*[100], maxlen=3)

    lr = 0.00001
    for iter in range(2000):
        loss = -(num_head * torch.log(theta) + num_tail * torch.log(1 - theta))
        loss.backward()
        with torch.no_grad():
            theta -= lr * theta.grad
            # print(f'{iter}: {theta}, {theta.grad}')
            recent.append(theta.grad.item())
            if all(map(lambda x: abs(x) < 1, recent)):
                break
        theta.grad.zero_()
    return theta.item()


if __name__ == '__main__':
    data = gen_coins(0.6, n=200)

    num_head = (data.detach() == 1).sum().item()
    num_tail = (data.detach() == 0).sum().item()

    print(num_head, num_tail)
    print(train(num_head, num_tail))

一点需要说明的是，在迭代过程中，我们保存最后三个导数的值，当最新的三个导数都很小时就退出迭代。

1

if all(map(lambda x: abs(x) < 1, recent))

运行代码，能发现最大化对数似然函数能很稳定的找到 \(\theta\)。

现在大家对于伯努利MLE有了一定了解，接着，我们来思考一下最大化似然函数方法是否随着观察次数的增多能不断逼近真实的 \(\theta^\star\)呢？

MLE \(\theta\) 估计的收敛性

\(\theta^\star=0.7\) 的情况下，我们来这样做试验，第一次做 n=1生成观察数据 \(X_{1}\)，第二次做 n=2生成观察数据 \(X_{2}\) \[ X_1,X_2, X_3, ..., X_N \] 对于每个数据集 \(X_i\) 通过最大似然方法求得估计的 \(\hat\theta\) \[ \hat\theta_1=MLE(X_1), \hat\theta_2=MLE(X_2), ..., \hat\theta_N=MLE(X_N) \] 将这些 \(\hat\theta_i\) 画出来，可以看到，随着 \(n \to \infty\)，\(\hat\theta_i \to \theta^\star=0.7\)

换一个角度来看一下，我们将 \(\hat\theta_i\) 数列按照顺序，离散化后再归一化比例，如下图画出来，红色的柱代表了最新的值 \(\hat\theta\)。可以发现，初始时候，\(\hat\theta\) 在较远离 0.7 的地方出现，随着 n 的增大，出现的位置比较接近 0.7。

MLE \(\theta\) 估计的偏差和方差

我们已经知道 MLE 方法可以通过观察数据推测出最有可能的 \(\hat\theta\)，由于观察数据 \(X\) 是伯努利过程产生的，具有随机性，那么 \(\hat\theta\) 可以看成是 \(\theta^\star\) 的随机变量。我们已经通过上面的试验知道随着试验次数的增大，我们的估计会越来越逼近真实值，现在的问题是对于固定的n，\(\hat\theta\) 的方差是多少，它的均值是否是无偏的呢？

带着这样的疑问，我们现在做如下试验：

固定 n=10，重复做实验，画出随着次数增多 \(\hat\theta\) 的分布，见图中绿色部分。同样的，红色是 n=80 不断试验的分布变换。

看的出来，随着试验次数的增多 - \(\hat\theta_{10}, \hat\theta_{80}\) 都趋近于正态分布

• \(\hat\theta_{10}\) 的分散度比 $ _{80}$ 要大，即方差要大

• \(\hat\theta_{10}, \hat\theta_{80}\) 的均值都在 0.7

上期 从零构建统计随机变量生成器之离散基础篇，我们从零出发构建了基于伯努利的基础离散分布，这一期我们来详细介绍广泛使用的 Inverse Transform Method（逆变换采样方法）。

• 从零构建统计随机变量生成器之 离散基础篇
• 从零构建统计随机变量生成器之 用逆变换采样方法构建随机变量生成器
• 深入 LeetCode 470 拒绝采样，状态转移图求期望和一道经典统计求期望题目
• 从零构建统计随机变量生成器之 正态分布 Box-Muller方法

逆变换采样方法

Inverse Transform Method 是最基础常见的方法，可用于离散分布和连续分布。常见的分布一般都能通过此方法生成，只需要随机变量CDF的解析表达式。假设随机变量 \(X\)，其CDF为 \(F^{-1}\)，则 Inverse Transform Method 仅有两步

1. 通过生成 [0, 1] 之间的均匀随机数 \(u\)
2. 代入 \(F^{-1}\) 即产生满足\(X\)分布的实例 \(x = F^{-1}(u)\)

离散例子

我们先举一个离散分布来直观感受一下其工作机制。有如下PMF的离散类别分布，范围在 [1, 5]。 \[ P(X = 1)=\frac{1}{15} \]

\[ P(X = 2)=\frac{2}{15} \]

\[ P(X = 3)=\frac{1}{5} \]

\[ P(X = 4)=\frac{4}{15} \]

\[ P(X = 5)=\frac{1}{3} \]

转换成CDF为

\[ P(X \leq 1)=\frac{1}{15} \]

\[ P(X \leq 2)=\frac{1}{15}+\frac{2}{15}=\frac{1}{5} \]

\[ P(X \leq 3)=\frac{1}{15}+\frac{2}{15}+\frac{1}{5}=\frac{6}{15} \]

\[ P(X \leq 4)=\frac{1}{15}+\frac{2}{15}+\frac{1}{5}+\frac{4}{15}=\frac{2}{3} \]

\[ P(X \leq 5)=\frac{1}{15}+\frac{2}{15}+\frac{1}{5}+\frac{4}{15}+\frac{1}{3}=1 \]

画出对应的CDF图

那么Inverse Transformation Method 的第一步，随机生成 0-1 之间的数 u，可以直观的认为是在 y 轴上生成一个随机的点 u。注意到5段竖虚线对应了5个离散的取值，它们的长度和为1，并且每一段长度代表了每个值的权重。因此，通过在 y 轴上的均匀采样可以生成给定PMF的 x 的分布。

离散分布的逆变换采样方法用数学公式可以表述为：找到第一个 x，其CDF的范围包括了 u，即

\[ F^{-1}(u)=\min \{x: F(x) \geq u\} \]

扩展到连续分布

有了离散类别分布的直观感受，扩展到连续分布也就不难理解了。类似于微积分中将连续空间做离散切分，再通过极限的方法，连续光滑函数在 y 轴上可以切分成长度为 \(\Delta u\) 的线段，那么生成的 x 值就是其近似值。随着 $ _{u } $，最终 $ x=F^{-1}(u) $ 即为满足要求的分布。

指数分布（连续）

以最为常见的指数分布为例，我们来看看具体的步骤。

我们知道指数分布的PDF如下

PDF 图为

计算CDF为

CDF 图

可以求得逆函数为

\[ x=F^{-1}(u)=-\frac{1}{\lambda} \ln (1-u) \]

由于 1-u 在 [0, 1] 范围上的随机数等价于 u，因此，x 的生成公式等价于

\[ x=-\frac{1}{\lambda} \ln (u) \]

实现代码

对应代码很简单

1
2
3
4
5
6

import random
from math import log2 as ln

def exp_gen(lambbda: float) -> float:
    u = random.random()
    return -ln(u) / lambbda

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/coutinous_exp_inv.py

类别分布（离散）

我们再来看基于类别分布 Inverse Transformation Method的其他离散分布例子。在从零构建统计随机变量生成器之离散基础篇中，我们已经介绍了类别分布（Categorical Distribution）的逆变换采样算法，同时还介绍了通过模拟 Bernoulli 实验来生成二项，几何，超几何分布的方法。在这一篇中，我们通过逆变换采样算法再来生成这些分布。

先回顾一下类别分布的逆变换采样实现。

给定如下的类别分布， $p = [0.2, 0.3, 0.1, 0.4] $

实现代码

类别分布的逆变换采样实现需要找到第一个大于 u 的元素的索引序号，在我们的实现中，将 $p = [0.2, 0.3, 0.1, 0.4] $ 转换成累计概率 $c = [0.2, 0.5, 0.6, 1] $ 后，由于 \(\vec c\) 数组是非递减的，因此我们可以用二分法代替线性查找，将时间复杂度降到 \(O(log(n))\)。下面的实现中直接调用 python bisect 函数即可。

1
2
3
4
5
6
7
8
9
10
11
12
13

import bisect
import random
from typing import List

def categorical(probs: List[float]) -> int:
    assert abs(sum(probs) - 1.0) < 0.001
    cum = probs.copy()

    for i in range(1, len(cum)):
        cum[i] = cum[i-1] + probs[i]

    u = random.random()
    return bisect.bisect(cum, u)

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_categorical.py

二项分布（离散）

二项分布（Binomial Distribution）有两个参数 n 和 p，表示伯努利实验做n次后成功的次数。图中为 n=6，p=0.5的二项分布。

概率质量函数（PMF）

实现代码

根据上面的PMF定义，我们将 [0, 6] 上的PMF计算出来，然后调用类别分布的逆变换采样实现即可：

1
2
3
4
5
6
7
8
9

from scipy.special import comb

from discrete_categorical import categorical
from math import pow


def binomial(n: int, p: float) -> int:
    pmf = [comb(n, k, exact=True) * pow(p, k) * pow(1-p, n-k) for k in range(0, n + 1)]
    return categorical(pmf)

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_binomial_inv.py

超几何分布（离散）

同样的，超几何分布（HyperGeometric Distribution）也可以如法炮制。

概率质量函数（PMF）

实现代码

1
2
3
4
5
6
7
8

from scipy.special import comb

from discrete_categorical import categorical

def hypergeometric(N: int, K_succ_num: int, n_trial_num: int) -> int:
    pmf = [comb(K_succ_num, k, exact=True) * comb(N - K_succ_num, n_trial_num - k, exact=True) / comb(N, n_trial_num, exact=True)
           for k in range(max(0, n_trial_num - (N - K_succ_num)), min(K_succ_num, n_trial_num) + 1)]
    return categorical(pmf)

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_hypergeometric_inv.py

几何分布（离散）

几何分布（Geometric Distribution）和上面的二项分布以及超几何分布不同的是，它的 support 是所有非负整数，因此，我们无法穷举计算所有 x 的概率。但是，我们可以通过将CDF 推出 Inverse CDF的解析表达式来直接实现。

概率质量函数（PMF）

\[ \operatorname{P}_\text{Geometric}(X=k)=(1-p)^{k-1} p \]

CDF

\[ F_X(x) = 1- (1-p)^x \]

Inverse CDF

反函数求得为 \[ F^{-1}(u) = \lfloor { log(1-u) \over log(1-p) }\rfloor \]

实现代码

1
2
3
4
5
6
7

import random
from math import floor
from math import log2 as ln

def geometric(p: float) -> int:
    u = random.random()
    return floor(ln(u) / ln(1-p))

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_geometric_inv.py

标准正态分布

一般，标准正态分布用Box-Muller 方法来生成，这个后续将做详细介绍。这里我们用 Schmeiser 提出的基于Inverse Transformation Method的近似方法来生成：

\[ X=F^{-1}(u) \approx \frac{u^{0.135}-(1-u)^{0.135}}{0.1975} \]

实现代码

1
2
3
4
5
6

import random

def normal():
    import math
    u = random.random()
    return (math.pow(u, 0.135) - math.pow(1-u, 0.135)) / 0.1975

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/coutinous_normal_apprx.py

http://www.columbia.edu/~ks20/4404-Sigman/4404-Notes-ITM.pdf

https://www.win.tue.nl/~marko/2WB05/lecture8.pdf

泊松分布

1
2
3
4
5
6
7
8
9
10
11
12
13

import random
from math import exp


def poisson(lambdda: float) -> int:
    total = 1.0
    i = 0
    threshold = exp(-1 * lambdda)
    while total >= threshold:
        u = random.random()
        total *= u
        i += 1
    return i - 1

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_poisson_inv.py

1
2
3
4
5
6
7
8
9
10

from numpy.random import exponential

def poisson(lambdda: float) -> int:
    total = 0.0
    i = 0
    while total <= lambdda:
        y = exponential(1)
        total += y
        i += 1
    return i - 1

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_poisson_from_exp.py

泊松分布

$ E_{1}, E_{2}, E_{3}, (1) $

\[ \mathbb{P}(K \geqslant k)=\mathbb{P}\left(E_{1}+\cdots+E_{k} \leqslant \lambda\right) \]

1
2
3
4
5
6
7
8
9
10

from numpy.random import exponential

def poisson(lambdda: float) -> int:
    total = 0.0
    i = 0
    while total <= lambdda:
        y = exponential(1)
        total += y
        i += 1
    return i - 1

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_poisson_from_exp.py

在本系列中，我们会从第一性原理出发，从零开始构建统计学中的常见分布的随机变量生成器，包括二项分布，泊松分布，高斯分布等。在实现这些基础常见分布的过程中，会展示如何使用统计模拟的通用技术，包括 inverse CDF，Box-Muller，分布转换等。本期通过伯努利试验串联起来基础离散分布并通过代码来实现这些分布的生成函数，从零开始构建的原则是随机变量生成器实现只依赖 random() 产生 [0, 1.0] 之间的浮点数，不依赖于其他第三方API来完成。

• 从零构建统计随机变量生成器之 离散基础篇
• 从零构建统计随机变量生成器之 用逆变换采样方法构建随机变量生成器
• 深入 LeetCode 470 拒绝采样，状态转移图求期望和一道经典统计求期望题目
• 从零构建统计随机变量生成器之 正态分布 Box-Muller方法

均匀分布（离散）

离散均匀分布（Discrete Uniform Distribution）的随机变量是最为基本的，图中为 [0, 6] 七个整数的离散均匀分布。算法实现为，使用 [0, 1] 之间的随机数 u，再将 u 等比例扩展到指定的整数上下界。

实现代码

1
2
3
4
5
6
7

import random
from math import floor

def uniform(a: int, b: int) -> int:
    assert a <= b
    u = random.random()
    return a + floor((b - a + 1) * u)

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_uniform.py

伯努利分布

伯努利分布（Bernoulli Distribution）是support为0或者1的离散分布，0和1可以看成失败和成功两种可能。伯努利分布指定了成功的概率p，例如，下图是 p=0.4 的伯努利分布。

伯努利分布随机数实现也很直接，将随机值 u 根据 p 决定成功或者失败。

实现代码

1
2
3
4
5
6

import random

def bernoulli(p: float) -> int:
    assert 0 <= p <= 1
    u = random.random()
    return 1 if u <= p else 0

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_bernoulli.py

类别分布

类别分布（Categorical Distribution）是在伯努利分布的基础上扩展到了多个点，每个点同样由参数指定了其概率，因此，其参数从 p 扩展到了向量 \(\vec p\)，如图所示为 $p = [0.2, 0.3, 0.1, 0.4] $ 时的类别分布。

实现代码

类别分布生成函数也扩展了伯努利分布的实现算法，将随机数 u 和累计概率向量作比较。在这个例子中， $p = [0.2, 0.3, 0.1, 0.4] $ 转换成 $c = [0.2, 0.5, 0.6, 1] $，再将 u 和 \(\vec c\)数组匹配，返回结果为第一个大于 u 的元素 index。实现上，我们可以以线性复杂度遍历数组，更好一点的方法是，用 python bisect函数通过二分法找到index，将时间复杂度降到 \(O(log(n))\)。

1
2
3
4
5
6
7
8
9
10
11
12
13

import bisect
import random
from typing import List

def categorical(probs: List[float]) -> int:
    assert abs(sum(probs) - 1.0) < 0.001
    cum = probs.copy()

    for i in range(1, len(cum)):
        cum[i] = cum[i-1] + probs[i]

    u = random.random()
    return bisect.bisect(cum, u)

Github 代码地址： https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_categorical.py

二项分布

二项分布（Binomial Distribution）有两个参数 n 和 p，表示伯努利实验做n次后成功的次数。图中为 n=6，p=0.5的二项分布。

实现代码

二项分布生成算法可以通过伯努利试验的故事来实现，即调用 n 次伯努利分布生成函数，返回总的成功次数。

1
2

def binomial(n: int, p: float) -> int:
    return sum(bernoulli(p) for _ in range(n))

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_binomial.py

概率质量函数（PMF）

几何分布

几何分布（Geometric Distribution）和伯努利实验的关系是：几何分布是反复伯努利实验直至第一次成功时的失败次数。如图，当成功概率 p=0.4时的几何分布。

实现代码

1
2
3
4
5
6
7

from discrete_bernoulli import bernoulli

def geometric(p: float) -> int:
    fail_num = 0
    while not bernoulli(p):
        fail_num += 1
    return fail_num

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_geometric.py

概率质量函数（PMF）

\[ \operatorname{Pr}_\text{Geometric}(X=k)=(1-p)^{k-1} p \]

负二项分布

负二项分布（Negative Binomial Distribution）是尝试伯努利试验直至成功 r 次的失败次数。

实现代码

1
2
3
4
5
6
7
8
9
10
11

from discrete_bernoulli import bernoulli

def negative_binomial(r: int, p: float) -> int:
    failures = 0
    while r:
        success = bernoulli(p)
        if success:
            r -= 1
        else:
            failures += 1
    return failures

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_nagative_binomial.py

概率质量函数（PMF）

超几何分布

超几何分布（HyperGeometric Distribution）的意义是从总数为N的集合抽取n次后成功的次数。具体来说，集合由K个表示成功的元素和N-K个表示失败的元素组成，并且抽取时没有替换（without replacement）情况下的成功次数。注意，超几何分布和二项分布的区别仅在于有无替换。

实现代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

from discrete_bernoulli import bernoulli

def hypergeometric(N: int, K_succ_num: int, n_trial_num: int) -> int:
    x = N - K_succ_num
    n_hit = 0
    while n_trial_num:
        hit = bernoulli(K_succ_num / (K_succ_num + x))
        n_hit += hit
        if hit:
            K_succ_num -= 1
        else:
            x -= 1
        if K_succ_num == 0:
            return n_hit
        n_trial_num -= 1
    return n_hit

Github 代码地址：

https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_hypergeometric.py

概率质量函数（PMF）

负超几何分布

负超几何分布（Negative Hypergeometric Distribution）的意义是从总数为N的集合中，无替换下抽取直至 r 次失败时，成功的次数。

实现代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

from discrete_bernoulli import bernoulli

def negative_hypergeometric(N: int, K_success_num: int, r_fail_times: int) -> int:
    fail_num = N - K_success_num
    succ_trials = 0
    while r_fail_times:
        success = bernoulli(K_success_num / (K_success_num + fail_num))
        if success:
            K_success_num -= 1
            succ_trials += 1
            if K_success_num == 0: # no more success elements
                return succ_trials
        else:
            fail_num -= 1
            r_fail_times -= 1
    return succ_trials

Github 代码地址： https://github.com/MyEncyclopedia/stats_simulation/blob/main/distrib_sim/discrete_negative_hypergeometric.py

概率质量函数（PMF）

伯努利试验总结

下表总结了上面四种和伯努利试验有关的离散分布的具体区别。