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

Question

假设我知道“成功”的概率是 P。我运行了 N 次测试，我看到 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