计算巨大数字的二项概率

Question

我想在Python中计算二项概率。我试着用公式来计算：

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

但我得到的一些概率是无限大的。我检查了一些情况下的值，其中p=无穷大。比如说，当n=450,000和k=17时，这个值必须大于1e302，而这是浮点数能处理的最大值。

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

这个方法会抽取numberOfTrials个样本，并计算valueOfInterest这个值出现的平均次数。

这样做不会出现无限大的值。但是，这样做是否合理呢？为什么用这种方法不会出现无限大，而直接计算概率时却会出现呢？

Answer 1

为了避免像零乘以无穷大这样的复杂情况，可以采用逐步相乘的方法，像这样。

def Pbinom(N,p,k):
    q=1-p
    lt1=[q]*(N-k)
    gt1=list(map(lambda x: p*(N-k+x)/x, range(1,k+1)))
    Pb=1.0
    while (len(lt1) + len(gt1)) > 0:
        if Pb>1:
            if len(lt1)>0:
                Pb*=lt1.pop()
            else:
                if len(gt1)>0:
                    Pb*=gt1.pop()
        else:
            if len(gt1)>0:
                Pb*=gt1.pop()
            else:
                if len(lt1)>0:
                    Pb*=lt1.pop()
    return Pb

Answer 2

我觉得你应该使用对数来进行所有的计算：

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 3

因为你在使用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 4

在对数域中进行计算，以便处理组合和指数函数，然后再将结果进行指数运算。

大致是这样的：

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 时会出现下溢的问题。不过，你仍然可以继续使用对数概率进行计算。