算法导论 练习5.4-4

一次聚会需要邀请多少人，才能让其中3人的生日很可能相同？

为了便于叙述，对题目作一个转化：

考虑将n粒豆子随机放入365个桶中，至少有一个桶内豆子数量不少于3的概率，其中豆子、桶均有区别

因为豆子和桶均有区别，可得全事件的数量为365^n

豆子放入完成后，桶内至多有一粒豆子的情况相当于从365个桶中选取n个，并将n粒不同豆子依次放入，因此事件数为：

仅有一个桶内有两粒豆子，其余均为一粒或没有。这里我们首先考虑从n粒豆子中选出两粒绑定在一起，进而转化为从365个桶中选取n-1个，将豆子依次放入，事件数为：

有且仅有i个桶内有两粒豆子，其余均为一粒或没有。这里有一点需要注意的是绑定豆子的顺序问题，即先绑定豆子[1, 2]，后绑定豆子[3, 4]与先绑定[3, 4]后绑定[1, 2]是相同的，因此在绑定后需除以排列数。事件数为：

从全事件中去除上述事件后，剩余事件占全事件的比例即为至少一桶内有不少于3粒豆子的概率，实现如下：

from math import factorial
import numpy as np

class Birthday:
    def C(self, a, b):
        return factorial(a) / factorial(a - b) / factorial(b)
    def A(self, a, b):
        return factorial(a) / factorial(a - b)
    def slove(self, n):
        if n < 3:
            return 0
        elif n > 730:
            return 1
        elif n < 366:
            pair = int(n / 2)
            p = (self.C(365, n) * self.A(n, n) + sum(self.C(365, n - i) * np.prod([self.C(n - 2 * k, 2) for k in range(i)]) / self.A(i, i) * self.A(n -i, n- i) for i in range(1, pair + 1))) / 365**n
            return 1 - p
        else:
            # 估计用不上，计算量太大
            pair = int(n / 2)
            return 1 - sum(self.C(365, n - i) * np.prod([self.C(n - 2 * k, 2) for k in range(i)]) / self.A(i, i) * self.A(n -i, n- i) for i in range(1, pair + 1)) / 365**n

if __name__ == '__main__':
    bir = Birthday()
    bir.slove(30)   # 0.0285
    bir.slove(50)   # 0.1264
    bir.slove(80)   # 0.4182
    bir.slove(87)   # 0.4995
    bir.slove(88)   # 0.5111
    bir.slove(120)  # 0.8280

可以看出，88人时，有三人生日相同的概率即超过0.5

此处答案与使用期望计算略有不同，个人认为是因为期望计算时包含了部分重复导致