Wie kann ich effizient berechnet die binomiale kumulative Verteilungsfunktion?

Sagen wir mal so, ich weiß, die Wahrscheinlichkeit eines "Erfolg" ist P. ich den test ausführen, N mal, und ich sehe S Erfolge. Der test ist vergleichbar mit dem werfen einer ungleich gewichteten Münze (vielleicht Köpfe ist ein Erfolg, tails ist ein Fehler).

Ich möchte wissen, die Ungefähre Wahrscheinlichkeit des Sehens entweder S Erfolge, oder eine Reihe von erfolgen weniger wahrscheinlich als S-Erfolge.

So zum Beispiel, wenn P ist 0,3, N-100 und ich bekommen, 20 Erfolge, ich bin auf der Suche nach der Wahrscheinlichkeit 20 oder weniger Erfolge.

Wenn auf der anderen hatte, P ist 0,3, N 100, und ich bekomme 40 Erfolge, ich bin auf der Suche nach der Wahrscheinlichkeit, 40 unsere weitere Erfolge.

Ich bin mir bewusst, dass dieses problem bezieht sich auf das finden der Fläche unter einer binomialen Kurve, jedoch:

Meine Mathe-fu ist nicht bis zu der Aufgabe der übersetzung dieser Kenntnisse in effizienten code
Ich verstehe zwar einer binomialen Kurve geben würde, ein exaktes Ergebnis, ich habe den Eindruck, es wäre grundsätzlich ineffizient. Eine schnelle Methode zum berechnen einer ungefähren Ergebnis, würde genügen.

Ich sollte betonen, dass diese Berechnung muss schnell sein, und sollte idealerweise definierbar mit standard-64-oder 128-bit-floating-point-Berechnung.

Ich bin auf der Suche nach einer Funktion, P, S, und N - und gibt eine Wahrscheinlichkeit. Als ich bin mehr vertraut mit dem code, als die mathematische Schreibweise ist, würde ich es vorziehen, dass alle Antworten beschäftigen pseudo-code oder code.

InformationsquelleAutor der Frage sanity | 2009-07-08

Exakten Binomialen Verteilung

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

Normalen Schätzen, gut für große 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

Poisson-Schätzung: für große n und kleine 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

InformationsquelleAutor der Antwort Unknown

6

War ich an einem Projekt, wo wir benötigt, um in der Lage sein, um die Berechnung der binomial CDF-in einer Umgebung, die nicht über eine Faktoren-oder gamma-Funktion definiert. Es dauerte ein paar Wochen, aber ich landete kommenden up mit dem folgenden Algorithmus berechnet die CDF-genau (d.h. keine Angleichung erforderlich). Python ist im Grunde so gut wie pseudocode, richtig?
```
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
```
Leistung skaliert mit x. Für kleine Werte von x, diese Lösung ist etwa eine Größenordnung schneller als scipy.stats.binom.cdf mit ähnlicher Leistung bei etwa x=10,000.

Ich gehe nicht in eine vollständige Herleitung dieses Algorithmus, da stackoverflow nicht unterstützt MathJax, aber der Schub ist es zunächst der Identifizierung der folgenden äquivalenz:
- Für alle k > 0, sp.misc.comb(n,k) == np.prod([(n-k+1)/k for k in range(1,k+1)])
Können wir umschreiben als:
- sp.misc.comb(n,k) == sp.misc.comb(n,k-1) * (n-k+1)/k
oder in eine log-space:
- np.log( sp.misc.comb(n,k) ) == np.log(sp.misc.comb(n,k-1)) + np.log(n-k+1) - np.log(k)
Weil die CDF ist eine Summierung von PMFs, wir können mit dieser Formulierung zur Berechnung der binomial-Koeffizienten (log ist b in der Funktion oben) für PMF_{x=i} aus den Koeffizienten berechnet wir für PMF_{x=i-1}. Dies bedeutet, dass wir tun können, alles in einer einzigen Schleife mit Akkumulatoren, und wir brauchen nicht zu berechnen, beliebige faktorielle!

Der Grund, die meisten der Berechnungen sind im Protokoll-Raum ist die Verbesserung der numerischen Stabilität des Polynoms, D. H. p^x und (1-p)^(1-x) haben das Potenzial zu werden, extrem groß oder extrem klein, was kann die Ursache für Rechenfehler.

EDIT: Ist das ein neuartiger Algorithmus? Ich habe herumgestöbert und ausschalten da, bevor ich dies geschrieben, und ich bin immer gefragt, ob ich schreiben soll, das ist mehr formal und senden Sie an eine Zeitschrift.

InformationsquelleAutor der Antwort David Marx
4

Ich denke, Sie wollen zu bewerten, die unvollständige beta-Funktion.

Gibt es eine schöne Umsetzung mit einer Fortsetzung Bruchteil Darstellung in "Numerical Recipes In C", 6. Kapitel: 'Spezielle Funktionen'.

InformationsquelleAutor der Antwort duffymo
4

Kann ich nicht ganz bürgen für die Effizienz, aber Scipy hat eine Modul für das
```
from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)
```
InformationsquelleAutor der Antwort mcstrother
2

Einen effizienten und vor allem numerisch stabilen Algorithmus ist in der Domäne der Bezier-Kurven verwendet in Computer-Aided Design. Es heißt de Casteljau-Algorithmus bewerten, die Bernstein-Polynome verwendet, um zu definieren, Bezier-Kurven.

Ich glaube, ich darf nur ein link pro Antwort so starten Sie mit Wikipedia - Bernstein-Polynome

Beachten Sie die sehr enge Beziehung zwischen der Binomialverteilung und der Bernstein-Polynome. Dann klicken Sie sich durch den link auf de Casteljau-Algorithmus.

Können sagen, ich weiß, die Wahrscheinlichkeit des Werfens einer Köpfe mit einer bestimmten Münze ist in P.
Was ist die Wahrscheinlichkeit, mich zu werfen
die Münze T mal und immer mindestens
S Köpfe?
- Setzen n = T
- Set beta[i] = 0 für i = 0, ... N - 1
- Set beta[i] = 1 für i = S, ... T
- Menge t = p
- Bewerten B(t) mit de Casteljau
oder an den meisten S-Köpfe?
- Setzen n = T
- Set beta[i] = 1 für i = 0, ... S
- Set beta[i] = 0 für i = S + 1, ... T
- Menge t = p
- Bewerten B(t) mit de Casteljau
Open-source-code wohl bereits existiert. NURBS-Kurven (Non-Uniform Rational B-spline-Kurven) sind eine Verallgemeinerung von Bézier-Kurven und sind weit verbreitet in CAD. Versuchen openNurbs (die Lizenz ist sehr liberal) oder nicht, die Open CASCADE (eine etwas weniger liberale und undurchsichtig-Lizenz). Beide toolkits sind in C++, obwohl, wenn ich mich Recht erinnere, .NET-Bindungen existieren.

InformationsquelleAutor der Antwort Paul Delhanty

Wenn Sie mit Python, keine Notwendigkeit, code-it-yourself. Scipy haben Sie abgedeckt:

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

InformationsquelleAutor der Antwort volodymyr

1

Aus dem Teil deiner Frage "immer mindestens N-Köpfe" Sie wollen die kumulative Binomialverteilung Funktion. Sehen http://en.wikipedia.org/wiki/Binomial_distribution für die Gleichung, die beschrieben wird als in Bezug auf die "regulierte unvollständige beta-Funktion" (wie schon beantwortet). Wenn Sie nur wollen, um zu berechnen, die Antwort, die zimmerreserviereung, ohne das Sie implementieren die gesamte Lösung selbst, die GNU Scientific Library liefert die Funktion: gsl_cdf_binomial_P und gsl_cdf_binomial_Q.

InformationsquelleAutor der Antwort
1

Den DCDFLIB Projekt hat C# - Funktionen (Wrapper für C-code) zu bewerten, viele CDF-Funktionen, einschließlich der binomial-Verteilung. Sie finden die original-C-und FORTRAN-code hier. Dieser code ist getestet und korrekt.

Wenn Sie wollen, um Ihren eigenen code schreiben, um nicht abhängig von einer externen Bibliothek, Sie könnte verwenden Sie die normale Annäherung an die Binomialverteilung erwähnt in anderen Antworten. Hier sind einige Hinweise auf wie gut die Näherung ist unter verschiedenen Umständen. Wenn Sie diesen Weg gehen und müssen code zum berechnen der normal-CDF, hier ist Python-code tun. Es ist nur etwa ein Dutzend Zeilen code, und könnte leicht portiert werden, um eine andere Sprache. Aber wenn Sie wollen, hohe Genauigkeit und effizienten code, sind Sie besser dran, unter Verwendung von Drittanbieter-code, wie DCDFLIB. Mehrere Mann-Jahre ging in die Produktion, die Bibliothek.

InformationsquelleAutor der Antwort John D. Cook
0

Versuchen diese eineverwendet in der GMP. Ein weiterer Verweis diese.

InformationsquelleAutor der Antwort lhf

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.

InformationsquelleAutor der Antwort ramakrishnareddy

Schreibe einen Kommentar

Du musst angemeldet sein, um einen Kommentar abzugeben.