18

「成功」の確率が 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 を取り、確率を返す関数を探しています。私は数学表記よりもコードに精通しているので、回答には疑似コードまたはコードを使用することをお勧めします。

4

10 に答える 10

32

正確な二項分布

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
于 2009-07-11T19:42:19.083 に答える
5

不完全なベータ関数を評価したいと思います。

「Numerical Recipes In C」の第 6 章「Special Functions」には、連分数表現を使用した優れた実装があります。

于 2009-07-08T01:38:53.053 に答える
4

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

from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)
于 2012-02-08T23:40:51.430 に答える
3

コンピュータ支援設計で使用されるベジエ曲線の領域には、効率的で、さらに重要なことに、数値的に安定したアルゴリズムが存在します。これは、ベジエ曲線を定義するために使用されるバーンスタイン多項式を評価するために使用されるde Casteljau のアルゴリズムと呼ばれます。

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

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

特定のコインで表が出る確率が P であることがわかっているとします。コインを T 回投げて、少なくとも S 表が出る確率は?

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

またはせいぜいS頭？

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

オープン ソース コードはおそらく既に存在します。 NURBS 曲線(不均一な有理 B スプライン曲線) は、ベジエ曲線を一般化したもので、CAD で広く使用されています。openNurbs (ライセンスは非常にリベラルです) を試すか、その Open CASCADE (ややリベラルで不透明なライセンス) に失敗します。どちらのツールキットも C++ ですが、IIRC と .NET バインディングが存在します。

于 2009-07-08T05:28:32.860 に答える
2

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
于 2016-04-08T02:07:04.900 に答える
1

あなたの質問の「少なくともS頭を得る」という部分から、累積二項分布関数が必要です。「正規化された不完全なベータ関数」の観点から説明されている方程式については、http://en.wikipedia.org/wiki/Binomial_distributionを参照してください (既に回答済み)。ソリューション全体を自分で実装する必要がなく、単に答えを計算したい場合は、GNU Scientific Library に関数 gsl_cdf_binomial_P および gsl_cdf_binomial_Q が用意されています。

于 2009-07-08T04:59:36.023 に答える
1

DCDFLIB プロジェクトには、二項分布を含む多くの CDF 関数を評価するためのC# 関数 (C コードのラッパー) があります。元の C および FORTRAN コードは、ここにあります。このコードは十分にテストされており、正確です。

外部ライブラリに依存しないように独自のコードを作成する場合は、他の回答で言及されている二項式の正規近似を使用できます。ここでは、さまざまな状況下で近似がどの程度優れているかについて、いくつかの注意事項を示します。そのルートに進み、通常の CDF を計算するためのコードが必要な場合は、それを行うためのPython コードを次に示します。これはわずか 12 行のコードであり、他の言語に簡単に移植できます。ただし、高精度で効率的なコードが必要な場合は、DCDFLIB などのサードパーティ コードを使用することをお勧めします。そのライブラリの作成には数年が費やされました。

于 2009-07-12T02:02:36.000 に答える
0

GMPで使用されているこれを試してください。別の参照はこれです。

于 2009-07-08T01:21:08.180 に答える