二項累積分布関数を効率的に計算するにはどうすればよいですか?

https://stackoverflow.com/questions/1095650

11-09-2019
|

質問

「成功」の確率が 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 の値が小さい場合、この解はこれよりも約 1 桁高速です。 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 の合計であるため、この公式を使用して二項係数を計算できます (その対数は b 上記の関数で)、PMF_{x=i-1} に対して計算した係数から PMF_{x=i} に対して計算します。これは、アキュムレータを使用して単一のループ内ですべてを行うことができ、階乗を計算する必要がないことを意味します。

ほとんどの計算が対数空間で行われる理由は、多項式項の数値安定性を向上させるためです。 p^x そして (1-p)^(1-x) 非常に大きいか非常に小さい可能性があり、計算エラーが発生する可能性があります。

編集：これは新しいアルゴリズムですか?私はこれを投稿する前から断続的に調べ続けてきましたが、これをもっと正式に書いてジャーナルに投稿すべきかどうか、ますます疑問に思っています。

私はあなたがhref="http://mathworld.wolfram.com/BinomialDistribution.html"が不完全ベータ関数をrel="noreferrer">

、第6章「Cにおける数値のレシピ」の連分数表現を使って素敵な実装があります：「特殊機能は、」

効率を完全に保証することはできませんが、Scipy にはこのためのモジュール

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

効率的で、さらに重要なことに、数値的に安定したアルゴリズムは、次の領域に存在します。 ベジェ曲線 コンピュータ支援設計で使用されます。いわゆる ド・カステルジョーのアルゴリズム を評価するために使用されます バーンスタイン多項式 ベジェ曲線を定義するために使用されます。

回答ごとに 1 つのリンクのみが許可されていると思うので、次から始めてくださいウィキペディア - バーンスタイン多項式

二項分布とバーンスタイン多項式の間には非常に密接な関係があることに注目してください。次に、de Casteljau のアルゴリズムに関するリンクをクリックします。

特定のコインで表が出る確率が P であることがわかっているとします。私がコインを投げて少なくともSの頭を手に入れる確率はどのくらいですか？

n = T を設定します
i = 0 の場合、beta[i] = 0 を設定します。S-1
i = S の場合、beta[i] = 1 を設定します。T
t = p に設定します
de Casteljau を使用して B(t) を評価する

それともせいぜいS頭くらいでしょうか？

n = T を設定します
i = 0 の場合、beta[i] = 1 を設定します。S
i = S + 1 の場合、beta[i] = 0 を設定します。T
t = p に設定します
de Casteljau を使用して B(t) を評価する

オープンソースコードはおそらくすでに存在します。 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

は、あなたの質問の一部からは、cummulative二項分布関数をしたい「とは、少なくともSヘッドを取得します」。である、式のために http://en.wikipedia.org/wiki/Binomial_distributionするを参照してください。（既に答えとして）「正則不完全ベータ関数」の用語であると記載。あなただけのソリューション全体を自分で実装することなく、答えを計算したい場合は、GNU科学ライブラリが機能を提供します。gsl_cdf_binomial_Pとgsl_cdf_binomial_Qを

DCDFLIBプロジェクトには、多くを評価するために、C＃の機能（Cコードのラッパー）を持っています二項分布を含むCDF関数。あなたは、元のCとFORTRANのコードここに見つけることができます。このコードは十分にテストし、正確である。

あなたは外部ライブラリに依存しないようにする独自のコードを書きたい場合は、

、あなたは他の回答で述べた二項に正規近似を使用することができます。ここにいくつかの注意事項がある近似は、様々な状況下でしているどのように良いです。あなたはそのルートを行くと、通常のCDFを計算するためのコードが必要な場合は、ここでそれを行うためののPythonコードするです。それだけで、コードのダースラインについてです、簡単に他の言語に移植することができます。あなたは、高精度かつ効率的なコードをしたい場合しかし、あなたはDCDFLIBのようなサードパーティのコードを使用して方がいいでしょう。いくつかの人年、そのライブラリを生産に入っています。

GMPに使用この1 して、試してみてください。別の参照は、こののです。

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.

ライセンス： CC-BY-SA と帰属

所属していません StackOverflow