假设我知道“成功”的概率是 P。我运行测试 N 次，我看到 S 成功。该测试类似于扔一枚重量不均匀的硬币（也许正面是成功的，反面是失败的）。
我想知道看到 S 次成功或成功次数少于 S 次成功的大致概率。
因此，例如，如果 P 为 0.3，N 为 100，并且我获得 20 次成功，我正在寻找获得 20 次或更少成功的概率。
另一方面，如果 P 为 0.3，N 为 100，并且我获得 40 次成功，我正在寻找获得 40 次更多成功的概率。
我知道这个问题与找到二项式曲线下的区域有关，但是：
我应该强调，这种计算必须很快，并且理想情况下应该可以通过标准的 64 位或 128 位浮点计算来确定。
我正在寻找一个接受 P、S 和 N 并返回概率的函数。由于我更熟悉代码而不是数学符号，因此我希望任何答案都使用伪代码或代码。
精确二项分布
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
我在一个项目中，我们需要能够在没有定义阶乘或伽马函数的环境中计算二项式 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，所以我不会详细介绍该算法，但它的重点是首先确定以下等价性：
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 的总和，所以我们可以使用这个公式根据b我们为 PMF_{x=i-1 计算的系数计算 PMF_{x=i} 的二项式系数（其对数在上面的函数中） }。这意味着我们可以使用累加器在单个循环内完成所有操作，并且我们不需要计算任何阶乘！
大多数计算在对数空间中进行的原因是为了提高多项式项的数值稳定性，即p^x有(1-p)^(1-x)可能非常大或非常小，这会导致计算错误。
编辑：这是一种新颖的算法吗？自从我发布这篇文章之前，我就一直在四处寻找，我越来越想知道我是否应该更正式地写这篇文章并将其提交给期刊。
我想你想评估不完整的 beta 函数。
在“C 中的数值配方”第 6 章：“特殊函数”中有一个很好的使用连分数表示的实现。
我不能完全保证效率，但 Scipy 有一个模块
from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)
在计算机辅助设计中使用的贝塞尔曲线领域中存在一种有效且更重要的是数值稳定的算法。它被称为de Casteljau 算法，用于评估用于定义贝塞尔曲线的伯恩斯坦多项式。
我相信每个答案只允许我一个链接，所以从Wikipedia - Bernstein Polynomials开始
请注意二项式分布和伯恩斯坦多项式之间非常密切的关系。然后点击进入 de Casteljau 算法的链接。
假设我知道用特定硬币投掷正面的概率是 P。我投掷硬币 T 次并获得至少 S 个正面的概率是多少？
或者最多S个头？
开源代码可能已经存在。 NURBS 曲线（Non-Uniform Rational B-spline Curves）是贝塞尔曲线的推广，在 CAD 中广泛使用。尝试使用 openNurbs（许可证非常自由）或不使用 Open CASCADE（不那么自由和不透明的许可证）。这两个工具包都在 C++ 中，但存在 IIRC、.NET 绑定。
如果您使用的是 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
从您的问题“至少获得 S 头”部分，您需要累积二项式分布函数。有关方程式，请参见http://en.wikipedia.org/wiki/Binomial_distribution，该方程式被描述为“正则化不完全 beta 函数”（如已回答）。如果您只想计算答案而不必自己实现整个解决方案，GNU 科学库提供了函数：gsl_cdf_binomial_P 和 gsl_cdf_binomial_Q。
DCDFLIB项目具有 C# 函数（围绕 C 代码的包装器）来评估许多 CDF 函数，包括二项分布。您可以在此处找到原始的 C 和 FORTRAN 代码。此代码经过良好测试且准确。
如果您想编写自己的代码以避免依赖外部库，则可以使用其他答案中提到的二项式的正常近似值。这里有一些关于在各种情况下近似值有多好的注释。如果你走那条路并且需要代码来计算正常的 CDF，这里是Python 代码。它只有十几行代码，可以很容易地移植到任何其他语言。但是如果你想要高精度和高效的代码，你最好使用像 DCDFLIB 这样的第三方代码。制作那个图书馆花了好几个人年。
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.