如何用Python模拟含阶乘的生日概率问题？代码遇RuntimeError求解

错误原因 #

原代码直接计算factorial(365)会生成远超浮点数表示范围的极大数值，触发RuntimeError（数值溢出）。

修复方案 #

方案1：用排列数函数替代阶乘计算 #

factorial(N)/factorial(N-k)本质是从365天中选k天的排列数，scipy提供的perm函数可直接计算该值，避免大阶乘溢出：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import perm

N = 365
k = np.arange(1, 41)

def fun(k):
    # perm(N, k) 直接计算排列数P(N,k) = N!/(N-k)!
    return 1 - perm(N, k) / (N ** k)

plt.plot(k, fun(k))
plt.xlabel('人数k')
plt.ylabel('至少两人同天生日的概率')
plt.title('生日概率问题')
plt.show()

方案2：逐步乘积计算（无额外依赖） #

如果不想依赖scipy，可通过累积乘积计算“所有人生日都不同”的概率，再用1减去该值，完全规避大数值运算：

import numpy as np
import matplotlib.pyplot as plt

N = 365
k = np.arange(1, 41)

# 计算所有人生日不同的概率：(365/365)*(364/365)*...*((365-k+1)/365)
prob_diff = np.cumprod((N - np.arange(k)) / N)
prob_same = 1 - prob_diff

plt.plot(k, prob_same)
plt.xlabel('人数k')
plt.ylabel('至少两人同天生日的概率')
plt.title('生日概率问题')
plt.show()

说明 #

两种方案都能正确输出结果，方案2更轻量化，仅依赖numpy。运行后会得到一条上升曲线，当k=23时概率超过50%，符合生日悖论的经典结论。

内容的提问来源于stack exchange，提问作者kang kang