有理概率的离散采样：进制展开与期望抽样次数

问题

每一轮可以投掷两颗 6 面骰子。以可能少的轮数拟合概率的事件？

这里投掷两颗骰子一共有种可能，相当于一个的等概率随机数，因此问题可以表述为：

只允许调用 random.randint(0, 35)，如何定义事件使得？

转换成代码描述就是：

def f():
    x = random.randint(0, 35)
    # 待实现
    return 1

cnt, N = 0, 1000
for _ in range(N):
    if f() == 1:
        cnt += 1

实现f使得当时，

最简单思路

如果我们手上不是36种可能，而是只有26种，那么我们随机到0-13的概率就是，随机到14-25的概率就是，因此我们只需要在投掷到26-35这另外10种情况时再投下一轮就可以了。

def f():
    while True:
        x = random.randint(0, 35)
        if x <= 13:
            return 1
        elif x <= 25:
            return 0
        # 投掷到26-35时，进入下一轮

但是，我们可以发现在这套方案里：

每轮结束（不需要再递归）的概率
（落在共 26 个格子）
需要继续递归的概率

令为调用 randint(0,35) 的轮数，则

这就是参数的几何分布（数“直到第一次成功所需的试次数”）。

几何分布的期望：

这版“落 26 个格子结束、落 10 个格子重掷”的方案，期望轮数是。

通用公式：

若每轮有个等可能结果，其中个结果会“重掷”，则

剩下的问题就是，怎么样减少所需要投掷的轮数？所谓最优方案，就是找到一个期望投掷轮数最少的方案。

最优方案：进制展开

我们研究问题的一般形式：使用随机数采样出目标概率，其中。

将用进制表示: 若的质因子均来自，即，则在 base- 下为有限小数; 若，则为纯循环小数；一般情形为混循环（有限前缀 + 循环尾部）。

也就是我们能在进制下把写成“小数”形式：

我们可以依次生成进制位（每次用一个均匀随机数），比较前缀，
这种方法的本质是 前缀比较法 (prefix comparison)：

随机数序列和概率展开的“数字串”逐位比较。
一旦随机串小于目标前缀，判定为 成功 (1)。
一旦大于，判定为 失败 (0)。
若完全相等，则继续生成下一位。

例子：与进制

目标是。
我们先把它写成进制：

计算过程如下：

乘以得，取整数部分。
取小数部分继续乘：，所以。
再取小数部分，。
以此类推…

所以的进制展开约为：

抽样算法

调用一次 rand = randint(0, 35)。
若 rand < d1，输出 1（事件发生）。
若 rand > d1，输出 0（事件不发生）。
若 rand == d1，则继续生成下一位并比较 d2，以此类推。

这个过程保证 精确等价 于概率。

代码实现

def f():
    p, q = 7, 13 # 目标概率
    n = 36       # base

    remainder = p
    while True:
        # 计算下一位 digit
        remainder *= n
        d, remainder = divmod(remainder, q)

        # 随机生成一位
        x = random.randint(0, n - 1)

        if x < d:
            return 1
        elif x > d:
            return 0
        # x == d 时继续循环，比较下一位 digit

期望抽样次数

对于无限小数部分，由于每一步都以的概率继续递归，期望抽样次数为一个收敛的几何级数。使用几何分布的期望：

对于，期望抽样次数约为。优于之前的。

总结

把分数写成 base- 展开，用逐位比较实现概率抽样。这样我们就能利用 等概率的面随机数 精确拟合任何有理概率。

若分母的质因子全部来自，展开是 有限小数；否则是 无限循环小数（可能有有限前缀）。
在逐位比较法中，事件在有限步内决定的概率趋近于 1，几乎必然有限终止。
对于 循环小数的情况，期望调用随机数的次数是一个收敛的几何级数，其极限为

$调用次数$

这是在 base- 下的 最优期望。
如果展开在位内终止（位有限小数），那么期望次数是