8

我想在 python 上计算二项式概率。我试图应用公式：

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

我得到的一些概率是无限的。我检查了一些 p=inf 的值。对于其中之一，n=450,000 和 k=17。该值必须大于 1e302，这是浮点数处理的最大值。

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

这将抽取 numberOfTrials 个样本并计算抽取 valueOfInterest 值的平均次数。

这不会提高任何无限价值。但是，这是一种有效的方法吗？为什么这种方式不会提高任何无限价值而计算概率呢？

4

4 回答 4

9

因为您使用的是 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])
于 2014-03-05T16:48:41.960 回答
7

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

像这样的东西：

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)

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

于 2014-03-05T15:43:09.930 回答
7

我认为您应该使用对数进行所有计算：

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))
于 2014-03-05T15:46:27.227 回答
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
于 2015-02-08T16:33:34.497 回答