计算大数的二项式概率

Question

我想计算python上的二项式概率。我试着运用这个公式：

probability = scipy.misc.comb(n,k)*(p**k)*((1-p)**(n-k))

我得到的一些可能性是无限的。我检查了一些p=inf的值，其中一个n=450000，k=17。此值必须大于1e302，后者是浮点处理的最大值。在

然后我尝试使用sum(np.random.binomial(n,p,numberOfTrials)==valueOfInterest)/numberOfTrials

这将绘制numberOfTrials样本，并计算利息值的平均绘制次数。在

这不会引发任何无穷大的值。然而，这是一种有效的方法吗？为什么这种方法不会产生任何无穷大的值，而计算概率呢？在

Answer 1

我认为你应该用对数来计算：

from scipy import special, exp, log
lgam = special.gammaln

def binomial(n, k, p):
    return exp(lgam(n+1) - lgam(n-k+1) - lgam(k+1) + k*log(p) + (n-k)*log(1.-p))

Answer 2

因为您使用的是scipy，我想我应该提到scipy已经实现了统计分布。还要注意，当n这么大时，二项分布很好地近似于正态分布（或泊松，如果p很小）。在

n = 450000
p = .5
k = np.array([17., 225000, 226000])

b = scipy.stats.binom(n, p)
print b.pmf(k)
# array([  0.00000000e+00,   1.18941527e-03,   1.39679862e-05])
n = scipy.stats.norm(n*p, np.sqrt(n*p*(1-p)))
print n.pdf(k)
# array([  0.00000000e+00,   1.18941608e-03,   1.39680605e-05])

print b.pmf(k) - n.pdf(k)
# array([  0.00000000e+00,  -8.10313274e-10,  -7.43085142e-11])

Answer 3

在对数域中计算组合函数和幂函数，然后将它们提升为指数。在

像这样：

combination_num = range(k+1, n+1)
combination_den = range(1, n-k+1)
combination_log = np.log(combination_num).sum() - np.log(combination_den).sum()
p_k_log = k * np.log(p)
neg_p_K_log = (n - k) * np.log(1 - p)
p_log = combination_log + p_k_log + neg_p_K_log
probability = np.exp(p_log)

由于数字太大而导致数值下溢/溢出。在使用n=450000和p = 0.5, k = 17的示例中，它返回p_log = -311728.4，即，最终概率的日志非常小，因此在使用np.exp时发生下溢。但是，您仍然可以使用对数概率。在