在学习了一些基本的统计变量生成法之后，这次我们来看看如何生成正态分布。它就是大名鼎鼎的 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 个输入随机数，输出一个随机采样。

在这篇文章中，我们从一道LeetCode 470 题目出发，通过系统地思考，引出拒绝采样（Reject Sampling）的概念，并探索比较三种拒绝采样地解法；接着借助状态转移图来定量计算采样效率；最后，我们利用同样的方法来解一道稍微复杂些的经典抛硬币求期望的统计面试题目。

Leetcode 470 用 Rand7() 实现 Rand10()

已有方法 rand7 可生成 1 到 7 范围内的均匀随机整数，试写一个方法 rand10 生成 1 到 10 范围内的均匀随机整数。

不要使用系统的 Math.random() 方法。

思考

• rand7()调用次数的 期望值 是多少 ?

• 你能否尽量少调用 rand7() ?

思维过程

我们已有 rand7() 等概率生成了 [1, 7] 中的数字，我们需要等概率生成 [1, 10] 范围内的数字。第一反应是调用一次rand7() 肯定是不够的，因为覆盖的范围不够。那么，就需要至少2次调用 rand7() 才能生成一次 rand10()，但是还要保证 [1, 10] 的数字生成概率相等，这个是难点。 现在我们先来考虑反问题，给定rand10() 生成 rand7()。这个应该很简单，调用一次 rand10() 得到 [1, 10]，如果是 8, 9, 10 ，则丢弃，重新开始，否则返回。想必大家都能想到这个朴素的方法，这种思想就是统计模拟中的拒绝采样（Reject Sampling）。

有了上面反问题的思考，我们可能会想到，rand7() 可以生成 rand5()，覆盖 [1, 5]的范围，如果将区间 [1, 10] 分成两个5个值的区间 [1, 5] 和 [6, 10]，那么 rand7() 可以通过先等概率选择区间 [1, 5] 或 [6, 10]，再通过rand7() 生成 rand5()就可以了。这个问题就等价于先用 rand7() 生成 rand2()，决定了 [1, 5] 还是 [6, 10]，再通过rand7() 生成 rand5() 。

解法一：rand2() + rand5()

我们来实现这种解法。下图为调用两次 rand7() 生成 rand10 数值的映射关系：横轴表示第一次调用，1，2，3决定选择区间 [1, 5] ，4，5，6选择区间 [6, 10]。灰色部分表示结果丢弃，重新开始 （注，若第一次得到7无需再次调用 rand7()）。

有了上图，我们很容易写出如下 AC 代码。

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

# AC
# Runtime: 408 ms, faster than 23.80% of Python3 online submissions for Implement Rand10() Using Rand7().
# Memory Usage: 16.7 MB, less than 90.76% of Python3 online submissions for Implement Rand10() Using Rand7().
class Solution:
    def rand10(self):
        while True:
            a = rand7()
            if a <= 3:
                b = rand7()
                if b <= 5:
                    return b
            elif a <= 6:
                b = rand7()
                if b <= 5:
                    return b + 5

标准解法：rand49()

从提交的结果来看，第一种解法慢于多数解法。原因是我们的调用 rand7() 的采样效率比较低，第一次有 1/7 的概率结果丢弃，第二次有 2/7的概率被丢弃。

如何在第一种解法的基础上提高采样效率呢？直觉告诉我们一种做法是降低上述 7x7 表格中灰色格子的面积。此时，会想到我们通过两次 rand7() 已经构建出来 rand49()了，那么再生成 rand10() 也规约成基本问题了。

下图为 rand49() 和 rand10() 的数字对应关系。

实现代码比较简单。注意，while True 可以去掉，用递归来代替。

1
2
3
4
5
6
7
8
9
10
11

# AC
# Runtime: 376 ms, faster than 54.71% of Python3 online submissions for Implement Rand10() Using Rand7().
# Memory Usage: 16.9 MB, less than 38.54% of Python3 online submissions for Implement Rand10() Using Rand7().
class Solution:
    def rand10(self):
        while True:
            a, b = rand7(),  rand7()
            num = (a - 1) * 7 + b
            if num <= 40:
                return num % 10 + 1

更快的做法

上面的提交结果发现标准解法在运行时间上有了不少提高，处于中等位置。我们继续思考，看看能否再提高采样效率。

观察发现，rand49() 有 9/49 的概率，生成的值被丢弃，原因是 [41, 49] 只有 9 个数，不足10个。倘若此时能够将这种状态保持下去，那么只需再调用一次 rand7() 而不是从新开始情况下至少调用两次 rand7()， 就可以得到 rand10()了。也就是说，当 rand49() 生成了 [41, 49] 范围内的数的话等价于我们先调用了一次 rand9()，那么依样画葫芦，我们接着调用 rand7() 得到了 rand63()。63 分成了6个10个值的区间后，剩余 3 个数。此时，又等价于 rand3()，循环往复，调用了 rand7() 得到了 rand21()，最后若rand21() 不幸得到21，等价于 rand1()，此时似乎我们走投无路，只能回到最初的状态，一切从头再来了。

改进算法代码如下。注意这次击败了 92.7%的提交。

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

# AC
# Runtime: 344 ms, faster than 92.72% of Python3 online submissions for Implement Rand10() Using Rand7().
# Memory Usage: 16.7 MB, less than 90.76% of Python3 online submissions for Implement Rand10() Using Rand7().
class Solution:
    def rand10(self):
        while True:
            a, b = rand7(),  rand7()
            num = (a - 1) * 7 + b
            if num <= 40: return num % 10 + 1
            a = num - 40
            b = rand7()
            num = (a - 1) * 7 + b
            if num <= 60: return num % 10 + 1
            a = num - 60
            b = rand7()
            num = (a - 1) * 7 + b
            if num <= 20: return num % 10 + 1

采样效率计算

通过代码提交的结果和大致的分析，我们已经知道三个解法在采样效率依次变得更快。现在我们来定量计算这三个解法。

我们考虑生成一个 rand10() 平均需要调用多少次 rand7()，作为采样效率的标准。

一种思路是可以通过模拟方法，即将上述每个解法模拟多次，然后用总的 rand7() 调用次数除以 rand10() 的生成次数即可。下面以解法三为例写出代码

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

# The rand7() API is already defined for you.
rand7_c = 0
rand10_c = 0

def rand7():
    global rand7_c
    rand7_c += 1
    import random
    return random.randint(1, 7)
    
def rand10():
    global rand10_c
    rand10_c += 1
    while True:
        a, b = rand7(), rand7()
        num = (a - 1) * 7 + b
        if num <= 40: return num % 10 + 1
        a = num - 40   # [1, 9]
        b = rand7()
        num = (a - 1) * 7 + b  # [1, 63]
        if num <= 60: return num % 10 + 1
        a = num - 60  # [1, 3]
        b = rand7()
        num = (a - 1) * 7 + b  # [1, 21]
        if num <= 20: return num % 10 + 1

if __name__ == '__main__':
    while True:
        rand10()
        print(f'{rand10_c} {round(rand7_c/rand10_c, 2)}')

运行代码，发现解法三的采样效率稳定在 2.19。

采样效率精确计算

计算解法二

为了精确计算三个解法的采样效率，我们通过代码得到对应的状态转移图来帮助计算。

例如，解法一可以对应到下图：初始状态 Start 节点中的 +2 表示经过此节点会产生 2次 rand7() 的代价。从 Start 节点有 40/49 的概率到达被接受状态 AC，有 9/49 概率到达拒绝状态 REJ。REJ 需要从头开始，则用虚线表示重新回到 Start节点，也就是说 REJ 的代价等价于 Start。注意，从某个节点出发的所有边之和为1。

有了上述状态转移关系图，我们令初始状态的平均代价为 \(C_2\)，则可以写成递归表达式，因为其中 REJ 的代价就是 \(C_2\)，即

\[ C_2 = 2 + (\frac{40}{49}\cdot0 + \frac{9}{49} C_2) \]

解得 \(C_2\)

\[ C_2 = 2.45 \]

计算解法一

同样的，对于另外两种解法，虽然略微复杂，也可以用同样的方法求得。

解法一的状态转移图为

递归方程表达式为

\[ C_1 = 1+\frac{3}{7} \cdot (1+\frac{5}{7} \cdot 0 + \frac{2}{7} \cdot C_1) \cdot2+ \frac{1}{7} \cdot (C_1) \]

解得 \(C_1\)

\[ C_1 = \frac{91}{30} \approx 3.03 \]

计算解法三

最快的解法三状态转移图为

递归方程表达式为

\[ C_3 = 2+\frac{40}{49} \cdot 0 + \frac{9}{49} (1+\frac{60}{63} \cdot 0 + \frac{3}{63} \cdot (1+\frac{20}{21} \cdot 0 + \frac{1}{21} \cdot C_3)) \]

解得 \(C_3\) \[ C_3 = \frac{329}{150} \approx 2.193 \]

至此总结一下，三个解法的平均代价为 \[ C_1 \approx 3.03 > C_2 = 2.45 > C_3 \approx 2.193 \] 这些值和我们通过模拟得到的结果一致。

稍难些的经典概率求期望题目

至此，LeetCode 470 我们已经分析透彻。现在，我们已经可以很熟练的将此类拒绝采样的问题转变成有概率的状态转移图，再写成递推公式去求平均采样的代价（即期望）。这里，如果大家感兴趣的话不妨再来看一道略微深入的经典统计概率求期望的题目。

问题：给定一枚抛正反面概率一样的硬币，求连续抛硬币直到两个正面（正面记为H，两个正面HH）的平均次数。例如：HTTHH是一个连续次数为5的第一次出现HH的序列。

分析问题画出状态转移图：我们令初始状态下得到第一个HH的平均长度记为 S，那么下一次抛硬币有 1/2 机会是 T，此时状态等价于初始状态，另有 1/2 机会是 H，我们记这个状态下第一次遇见HH的平均长度为 H（下图蓝色节点）。从此蓝色节点出发，当下一枚硬币是H则结束，是T是返回初始状态。于是构建出下图。

这个问题稍微复杂的地方在于我们有两个未知状态互相依赖，但问题的本质和方法是一样的，分别从 S 和 H 出发考虑状态的概率转移，可以写成如下两个方程式：

解得

因此，平均下来，需要6次抛硬币才能得到 HH，这个是否和你直觉的猜测一致呢？

这个问题还可以有另外一问，可以作为思考题让大家来练习一下：第一次得到 HT 的平均次数是多少？这个是否和 HH 一样呢？

导读：极大似然估计(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）

伯努利试验总结

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

Leetcode 1227 是一道有意思的概率题，本篇将从多个角度来讨论这道题。题目如下

有 n 位乘客即将登机，飞机正好有 n 个座位。第一位乘客的票丢了，他随便选了一个座位坐下。 剩下的乘客将会： 如果他们自己的座位还空着，就坐到自己的座位上， 当他们自己的座位被占用时，随机选择其他座位，第 n 位乘客坐在自己的座位上的概率是多少？

示例 1： 输入：n = 1 输出：1.00000 解释：第一个人只会坐在自己的位置上。

示例 2： 输入: n = 2 输出: 0.50000 解释：在第一个人选好座位坐下后，第二个人坐在自己的座位上的概率是 0.5。

提示： 1 <= n <= 10^5

假设规模为n时答案为f(n)，一般来说，这种递推问题在数学形式上可能有关于n的简单数学表达式（closed form），或者肯定有f(n)关于f(n-k)的递推表达式。工程上，我们可以通过通过多次模拟即蒙特卡罗模拟来算得近似的数值解。

Monte Carlo 模拟发现规律

首先，我们先来看看如何高效的用代码来模拟。根据题意的描述过程，直接可以写出下面代码。seats为n大小的bool 数组，每个位置表示此位置是否已经被占据。然后依次给第i个人按题意分配座位。注意，每次参数随机数范围在[0,n-1]，因此，会出现已经被占据的情况，此时需要再次随机，直至分配到空位。

暴力直接模拟

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

def simulate_bruteforce(n: int) -> bool:
    """
    Simulates one round. Unbounded time complexity.
    :param n: total number of seats
    :return: True if last one has last seat, otherwise False
    """

    seats = [False for _ in range(n)]

    for i in range(n-1):
        if i == 0:  # first one, always random
            seats[random.randint(0, n - 1)] = True
        else:
            if not seats[i]:  # i-th has his seat
                seats[i] = True
            else:
                while True:
                    rnd = random.randint(0, n - 1) # random until no conflicts
                    if not seats[rnd]:
                        seats[rnd] = True
                        break
    return not seats[n-1]

运行上面的代码来模拟 n 从 2 到10 的情况，每种情况跑500次模拟，输出如下

1
2
3
4
5
6
7
8
9
10

1 => 1.0
2 => 0.55
3 => 0.54
4 => 0.486
5 => 0.488
6 => 0.498
7 => 0.526
8 => 0.504
9 => 0.482
10 => 0.494

发现当 n>=2 时，似乎概率都是0.5。

标准答案

其实，这道题的标准答案就是 n=1 为1，n>=2 为0.5。下面是 python 3 标准答案。本篇后面会从多个角度来探讨为什么是0.5 。

1
2
3

class Solution:
    def nthPersonGetsNthSeat(self, n: int) -> float:
        return 1.0 if n == 1 else 0.5

O(n) 改进算法

上面的暴力直接模拟版本有个最大的问题是当n很大时，随机分配座位会产生大量冲突，因此，最坏复杂度是没有任何上限的。解决方法是每次发生随机分配时保证不冲突，能直接选到空位。下面是一种最坏复杂度O(n)的模拟过程，seats数组初始话成 0，1，...，n-1，表示座位号。当第i个人登机时，seats[i:n] 的值为他可以选择的座位集合，而seats[0:i]为已经被占据的座位集合。由于[i: n]是连续空间，产生随机数就能保证不冲突。当第i个人选完座位时，将他选中的seats[k]和seats[i] 交换，保证第i+i个人面临的seats[i+1:n]依然为可选座位集合。

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

def simulate_online(n: int) -> bool:
    """
    Simulates one round of complexity O(N).
    :param n: total number of seats
    :return: True if last one has last seat, otherwise False
    """

    seats = [i for i in range(n)]

    def swap(i, j):
        tmp = seats[i]
        seats[i] = seats[j]
        seats[j] = tmp

    # for each person, the seats array idx available are [i, n-1]
    for i in range(n-1):
        if i == 0:  # first one, always random
            rnd = random.randint(0, n - 1)
            swap(rnd, 0)
        else:
            if seats[i] == i:  # i-th still has his seat
                pass
            else:
                rnd = random.randint(i, n - 1)  # selects idx from [i, n-1]
                swap(rnd, i)
    return seats[n-1] == n - 1

递推思维

这一节我们用数学递推思维来解释0.5的解。令f(n) 为第 n 位乘客坐在自己的座位上的概率，考察第一个人的情况（first step analysis），有三种可能

1. 第一个人选了第一个即自己的座位，那么最后一个人一定能保证坐在自己的座位。
2. 第一个人选了最后一个人的座位，无论中间什么过程，最后一个人无法坐到自己座位
3. 第一个人选了第i个座位，(1<i<n)，那么第i个人前面的除了第一个外的人都会坐在自己位置上，第i个人由于没有自己座位，随机在剩余的座位1，座位 [i+1,n] 中随机选择，此时，问题转变为f(n-i+1)，如下图所示。

通过上面分析，得到概率递推关系如下

\[ f(n) = \begin{align*} \left\lbrace \begin{array}{r@{}l} 1 & & p=\frac{1}{n} \quad \text{选了第一个位置} \\\\\\ f(n-i+1) & & p=\frac{1}{n} \quad \text{选了第i个位置，1<i<n} \\\\\\ 0 & & p=\frac{1}{n} \quad \text{选了第n个位置} \end{array} \right. \end{align*} \]

即f(n)的递推式为： \[ f(n) = \frac{1}{n} + \frac{1}{n} \times [ f(n-1) + f(n-2) + ...+ f(2)], \quad n>=2 \] 同理，f(n+1)递推式如下 \[ f(n+1) = \frac{1}{n+1} + \frac{1}{n+1} \times [ f(n) + f(n-1) + ...+ f(2)] \] \((n+1)f(n+1) - nf(n)\) 抵消 \(f(n-1) + ...f(2)\) 项，可得 \[ (n+1)f(n+1) - nf(n) = f(n) \] 即 \[ f(n+1) = f(n) = \frac{1}{2} \quad n>=2 \]

用数学归纳法也可以证明 n>=2 时 f(n)=0.5。

简化的思考方式

我们再仔细思考一下上面的第三种情况，就是第一个人坐了第i个座位，1<i<n，此时，程序继续，不产生结果，直至产生结局1或者2，也就是case 1和2是真正的结局节点，它们产生的概率相同，因此答案是1/2。

从调用图可以看出这种关系，由于中间节点 f(4)，f(3)，f(2)生成Case 1和2的概率一样，因此无论它们之间是什么关系，最后结果都是1/2.

知乎上有个很形象的类比理解方式

考虑一枚硬币，正面向上的概率为 1/n，反面也是，立起来的概率为 (n-2)/n 。我们规定硬币立起来重新抛，但重新抛时，n会至少减小1。求结果为反面的概率。这样很显然结果为 1/2 。

这里，正面向上对应Case 2，反面对应Case 1。

这种思想可以写出如下代码，seats为 n 大小的bool 数组，当第i个人（0<i<n）发现自己座位被占的话，此时必然seats[0]没有被占，同时seats[i+1:]都是空的。假设seats[0]被占的话，要么是第一个人占的，要么是第p个人（p<i）坐了，两种情况下乱序都已经恢复了，此时第i个座位一定是空的。

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

def simulate(n: int) -> bool:
    """
    Simulates one round of complexity O(N).
    :param n: total number of seats
    :return: True if last one has last seat, otherwise False
    """

    seats = [False for _ in range(n)]

    for i in range(n-1):
        if i == 0:  # first one, always random
            rnd = random.randint(0, n - 1)
            seats[rnd] = True
        else:
            if not seats[i]:  # i-th still has his seat
                seats[i] = True
            else:
                # 0 must not be available, now we have 0 and [i+1, n-1],
                rnd = random.randint(i, n - 1)
                if rnd == i:
                    seats[0] = True
                else:
                    seats[rnd] = True
    return not seats[n-1]