Site icon Thoughts on Computing

Analytical Combinatorics:吃饼干问题

知乎被邀请回答了一个问题,和之前探索过的类似。问题如下:

10个不同的盒子,每个盒子里装有10块饼干。每次随机选择一个盒子,然后吃掉其中的一块饼干。直到第一次出现了任何一个空盒子就结束。求吃的饼干数目的期望。

这种问题是求第一次遇见的步长问题。我们可以转换为下面两步来求:

  1. 求给定某个n, 概率 \text{Prob}_n(n\text{块的时候任务尚且没有完成})
  2. 然后求和就得到了答案\sum_{n=0}^{\infty} \text{Prob}_n (就是第一次出现步长的期望)

原理很简单:

将问题抽象一点,有N个盒子,每个盒子有M块饼干。

因此这个问题就转换成了,特定步长任务没有结束的概率。吃n块饼干,所有可能的序列是N^n,求概率的最后一步,就只剩下求分子了,也就是所有序列中,没有搞定的序列有多少。这个计数我们用Analytical Combinatorics。

吃饼干的序列是一个字符串,字符集合是盒子的编号:c_1c_2c_3\dots c_n是吃饼干序列,其中c_i \in [1,N]。这个序列的等价模型就是N个盒子,每个盒子里是序列的下标的集合,表明第i块饼干来自这个盒子。由于我们研究的是没有搞定,故而,不能吃光盒子里的M块饼干,因此用符号方法,写出这个组合结构是:

\text{SEQ}_N( \sum_{k=0}^{M-1} \text{SET}_k(\bullet) )

这样就可以写出它对应的生成函数:

\Big( \sum_{k=0}^{M-1} \frac{z^k}{k!} \Big)^N

然后这个问题的解就是:

\sum_{n=0}^{N(M-1)} \frac{1}{N^n} n! [z^n] \Big( \sum_{k=0}^{M-1} \frac{z^k}{k!} \Big)^N

用一些简单的情况验证:

写一段Python程序计算这个式子也仿真一下:

import sys
from math import factorial as fact
from sympy import symbols, series
from collections import defaultdict
from random import randint


def compute(M, N):
    x = symbols('x')
    max_n = N*(M-1) + 1
    s = (sum(x**i/fact(i) for i in range(M)))**N
    taylor_expansion = series(s, x, 0, max_n)
    r = 0
    for i in range(max_n):
        coeff = taylor_expansion.coeff('x', i)
        r += (fact(i) * coeff / (N**i))
    return r

def sim(M, N, loop=10000):
    return sum(sim_one_loop(M, N) for i in range(loop)) / loop

def sim_one_loop(M, N):
    dt = defaultdict(int)
    n = 0
    while True:
        i = randint(1, N)
        dt[i] += 1
        n += 1
        if dt[i] == M:
            return n

print(float(compute(10,10)))
print(sim(10,10, 5000000))

跑出来的结果如下,非常接近。

57.47138382169789
57.4659616

精准数字是:

17959807444280590025974644911503703443990782567687928945219550218930752367809149/312500000000000000000000000000000000000000000000000000000000000000000000000000

关于有没有进一步办法化简,我没有到。

参考资料:

Exit mobile version