如何用python模拟带阶乘的概率问题

Question

我想模拟一个概率问题，但是代码出现RunTimeError

房间里有k个人。假设每个 一个人的生日很可能是一年中 365 天中的任何一天（我们排除 2 月 29 日），并且人们的生日是独立的（我们假设没有 房间里的双胞胎）。小组中有两个或更多人的概率是多少 生日一样吗？

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

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

def fun(k):
    return 1 - factorial(N)/N**k/factorial(N-k)

plt.plot(k, fun(k))

我想成功模拟概率问题

Answer 1

首先，你的功能是错误的。请检查第 2.1 节。以下参考资料：

https://users.cs.utah.edu/~jeffp/teaching/cs5955/L2-StatisticalPrinciples.pdf

然后，相应地调整您的功能：

    def fun(k):
    return 1 - (1-1/N)**(k*(k-1)/2)

剧情：

    plt.plot(k, fun(k))

结果：

Answer 2

从在 SciPy 中计算生日悖论 我们被警告说，如果直接实现 N 和 k 的适度大值将导致公式溢出，从而导致函数 fun 运行时错误。诀窍是使用对数来避免上溢或下溢。

SciPy 没有对数阶乘函数，但有对数伽玛函数，所以我们用它代替。

修改后的代码

import numpy as np
import matplotlib.pyplot as plt

from numpy.lib.scimath import log
from numpy import exp
from scipy.special import gammaln

def fun(k):
  # Use log of gamma since there is n
  return 1 - exp( gammaln(N+1) - gammaln(N-k+1) - k*log(N) )
  # replaces 1 - factorial(N)/N**k/factorial(N-k)
  
N = 365
k = np.arange(1, 41)

# Display
plt.plot(k, fun(k))
plt.show()

结果