我怎样才能有效地计算二项式累积分布函数？(How can I efficiently calcula

比方说，我知道，“成功”的概率为P.运行测试了N次了，我见S成功。 该试验类似于投掷的不均匀加权硬币（也许头是成功的，尾部是失败的）。

我想知道，看到既无S的成功，或者一些成功除S成功的可能性较低的近似概率。

因此，举例来说，如果P是0.3，N是100，我也得到20次成功，我在寻找获得20次或更少的成功的概率。

如果，另一方面hadn，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类似的性能。

我不会进入这个算法的完整推导，因为计算器不支持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是保偏光纤的总和，我们可以使用该制剂来计算二项式系数（日志其中是b在上面的功能），用于PMF_ {X = 1}从我们对PMF_ {计算出的系数X =异1}。 这意味着我们可以用蓄电池做一个循环内的一切，我们不需要任何计算阶乘！

大多数计算都是在日志空间做的原因是为了提高多项式的项的数值稳定性，即p^x和(1-p)^(1-x)有潜力成为非常大或非常小，这可导致计算错误。

编辑：这是一种新的算法？ 我已经和关闭闲逛，因为我张贴在此之前，和我越来越想知道如果我应该写这件事更正式，并提交了一份期刊。

Answer 3:

我认为，要评估不完整的测试功能 。

有一个很好的实现使用“数字食谱在C”一个连分数表示，第6章：“特殊功能”。

Answer 4:

我不能完全担保的效率，但SciPy的具有该模块

from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)

Answer 5:

一个高效，更重要的是，数值稳定的算法存在于计算机辅助设计中使用Bezier曲线的领域。 这就是所谓的去使用应用Casteljau算法来评估用于定义Bezier曲线的Bernstein多项式 。

我相信，我只允许每一个答案链接，以便下手维基百科-伯恩斯坦多项式

注意二项分布和Bernstein多项式之间的关系非常密切。 然后通过点击上德卡斯特里奥算法的链接。

可以说，我知道扔头与特定硬币的概率是P.什么是我的抛硬币t次并获得至少头顶的概率是多少？

• n设置= T
• 设置的β[I] = 0，对于i = 0，...的S - 1
• 设置的β[I]，其中i = S = 1，...Ť
• 集T = P
• 评价使用去应用Casteljau B（t）的

或最多头顶？

• n设置= T
• 设置的β[I] = 1对于i = 0，内容S
• 设置的β[I] = 0对于i = S + 1，...Ť
• 集T = P
• 评价使用去应用Casteljau B（t）的

开放的源代码可能已经存在。 NURBS曲线 （非均匀有理B样条曲线）是Bezier曲线的概括和广泛应用于CAD。 尝试openNurbs（许可非常宽松），或做不到这一点打开CASCADE（有略显宽松，不透明的许可证）。 这两个工具包是在C ++中，虽然存在IIRC，.NET绑定。

Answer 6:

如果您在使用Python，无需自己编写的。 SciPy的一应俱全：

from scipy.stats import binom
# probability that you get 20 or less successes out of 100, when p=0.3
binom.cdf(20, 100, 0.3)
>>> 0.016462853241869434

# probability that you get exactly 20 successes out of 100, when p=0.3
binom.pmf(20, 100, 0.3)
>>> 0.0075756449257260777

Answer 7:

从你的问题的一部分“至少让头顶”你想要的累积性二项分布函数。 见http://en.wikipedia.org/wiki/Binomial_distribution的方程，其被描述为在“正则化不完全β函数”的术语是（如已经回答）。 如果你只是想算出答案，而不必自己实现整个解决方案中，GNU科学图书馆提供的功能：gsl_cdf_binomial_P和gsl_cdf_binomial_Q。

Answer 8:

该DCDFLIB项目有C＃函数（约C代码封装），以评估许多CDF功能，包括二项分布。 你可以找到原来的C和Fortran代码在这里 。 此代码是经过充分测试和准确。

如果你想编写自己的代码，以避免依赖于外部库，你可以使用正常的近似在其他的答案中提到的二项式。 下面是一些关于笔记逼近有多好各种情况下。 如果你走这条路线，需要代码来计算正常CDF，这里的Python代码做这件事。 它大约只有十几行代码，并可以很容易地移植到任何其他语言。 但是，如果你想要高精确度和高效率的代码，你就要去使用像DCDFLIB第三方的代码更好。 几个人一年走进生产该库。

Answer 9:

试试这一个 ，在GMP使用。 另一个引用是这个 。

Answer 10:

import numpy as np
np.random.seed(1)
x=np.random.binomial(20,0.6,10000) #20 flips of coin,probability of 
                                 heads percentage and 10000 times 
                                  done.
sum(x>12)/len(x)

The output is 41% of times we got 12 heads.

文章来源: How can I efficiently calculate the binomial cumulative distribution function?