如何有效地计算二项式累积分布函数？

Question

让我们说我知道"成功"的概率是P.我进行了N次测试,我看到了S的成功.测试类似于抛掷不均匀加权的硬币(也许头部是成功的,尾部是失败的).

我想知道看到S成功的大概概率,或者比S成功的可能性要小得多.

因此,例如,如果P为0.3,N为100,并且我获得20次成功,我正在寻找获得20次或更少次成功的概率.

如果,另一方面没有,P是0.3,N是100,我获得40次成功,我正在寻找获得40次成功的可能性.

我知道这个问题与找到二项式曲线下的区域有关,但是:

1. 我的数学不能完成将这些知识转化为有效代码的任务
2. 虽然我理解二项式曲线会给出精确的结果,但我得到的印象是它本身效率低下.计算近似结果的快速方法就足够了.

我应该强调这个计算必须很快,理想情况下应该可以通过标准的64或128位浮点计算来确定.

我正在寻找一个带P,S和N的函数,并返回一个概率.由于我比代码更熟悉数学符号,我更喜欢任何答案都使用伪代码或代码.

Answer 1

精确二项分布

def factorial(n): 
    if n < 2: return 1
    return reduce(lambda x, y: x*y, xrange(2, int(n)+1))

def prob(s, p, n):
    x = 1.0 - p

    a = n - s
    b = s + 1

    c = a + b - 1

    prob = 0.0

    for j in xrange(a, c + 1):
        prob += factorial(c) / (factorial(j)*factorial(c-j)) \
                * x**j * (1 - x)**(c-j)

    return prob

>>> prob(20, 0.3, 100)
0.016462853241869437

>>> 1-prob(40-1, 0.3, 100)
0.020988576003924564

正常估计,适合大n

import math
def erf(z):
        t = 1.0 / (1.0 + 0.5 * abs(z))
        # use Horner's method
        ans = 1 - t * math.exp( -z*z -  1.26551223 +
                                                t * ( 1.00002368 +
                                                t * ( 0.37409196 + 
                                                t * ( 0.09678418 + 
                                                t * (-0.18628806 + 
                                                t * ( 0.27886807 + 
                                                t * (-1.13520398 + 
                                                t * ( 1.48851587 + 
                                                t * (-0.82215223 + 
                                                t * ( 0.17087277))))))))))
        if z >= 0.0:
                return ans
        else:
                return -ans

def normal_estimate(s, p, n):
    u = n * p
    o = (u * (1-p)) ** 0.5

    return 0.5 * (1 + erf((s-u)/(o*2**0.5)))

>>> normal_estimate(20, 0.3, 100)
0.014548164531920815

>>> 1-normal_estimate(40-1, 0.3, 100)
0.024767304545069813

泊松估计:适用于大n和小p

import math

def poisson(s,p,n):
    L = n*p

    sum = 0
    for i in xrange(0, s+1):
        sum += L**i/factorial(i)

    return sum*math.e**(-L)

>>> poisson(20, 0.3, 100)
0.013411150012837811
>>> 1-poisson(40-1, 0.3, 100)
0.046253037645840323

Answer 2

我在一个项目中,我们需要能够在没有定义阶乘或伽玛函数的环境中计算二项式CDF.我花了几个星期,但我最终提出了以下算法,它准确地计算了CDF(即没有必要的近似值).Python基本上和伪代码一样好吧？

import numpy as np

def binomial_cdf(x,n,p):
    cdf = 0
    b = 0
    for k in range(x+1):
        if k > 0:
            b += + np.log(n-k+1) - np.log(k) 
        log_pmf_k = b + k * np.log(p) + (n-k) * np.log(1-p)
        cdf += np.exp(log_pmf_k)
    return cdf

性能随x而变化.对于小的x值,该解决方案比scipy.stats.binom.cdf在x = 10,000附近的类似性能快约一个数量级.

我不打算完全推导出这种算法,因为stackoverflow不支持MathJax,但它的主旨是首先确定以下等价:

• 对于所有k> 0, sp.misc.comb(n,k) == np.prod([(n-k+1)/k for k in range(1,k+1)])

我们可以改写为:

• sp.misc.comb(n,k) == sp.misc.comb(n,k-1) * (n-k+1)/k

或在日志空间中:

• np.log( sp.misc.comb(n,k) ) == np.log(sp.misc.comb(n,k-1)) + np.log(n-k+1) - np.log(k)

因为CDF是PMF的总和,我们可以使用这个公式来计算bPMF_ {x = i} 的二项式系数(其中的对数在上面的函数中)来自我们为PMF_ {x = i-1计算的系数}.这意味着我们可以使用累加器在单个循环内完成所有操作,而且我们不需要计算任何因子!

大部分的计算都是在日志空间做的原因是为了提高多项式项的数值稳定性,也就是p^x和(1-p)^(1-x)必须非常大或非常小,这可能会导致计算错误的可能性.

编辑:这是一个新颖的算法？在我发布这篇文章之前,我一直在四处寻找,我越来越想知道我是否应该更正式地写这篇文章并将其提交给期刊.

Answer 3

我想你想要评估不完整的beta函数.

在"NC中的数字配方"第6章:"特殊功能"中使用连续分数表示有一个很好的实现.