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

假设我知道“成功”的概率为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并返回概率的函数。因为我比数学符号更熟悉代码,所以我希望任何答案都使用伪代码或代码。

回答:

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

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

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

以上是 如何有效计算二项式累积分布函数? 的全部内容, 来源链接: utcz.com/qa/433230.html

回到顶部