Wie kann ich die binomische Verteilungsfunktion effizient berechnen?

StackOverflow https://stackoverflow.com/questions/1095650

Frage

Lassen Sie uns sagen, dass ich die Wahrscheinlichkeit eines „Erfolg“ wissen, ist P. betreibe ich den Test N-mal, und ich sehe S Erfolge. Der Test ist ähnlich eine ungleich gewichteten Münze zu werfen (vielleicht Köpfe ist ein Erfolg, Schwanz ist ein Fehler).

Ich möchte sehen entweder S Erfolge die ungefähre Wahrscheinlichkeit kennen, oder eine Reihe von Erfolgen weniger wahrscheinlich als S Erfolge.

So zum Beispiel, wenn P 0,3, N 100, und ich bekomme 20 Erfolge, ich suche für die Wahrscheinlichkeit des Erhaltens 20 oder weniger Erfolge.

Wenn auf der anderen hadn, P 0,3, N 100, und ich bekomme 40 Erfolge, ich suche für die Wahrscheinlichkeit des Erhaltens 40 unserer mehr Erfolge.

Ich bin mir bewusst, dass dieses Problem betrifft den Bereich unter einer binomischen Kurve zu finden, aber:

  1. Meine Mathe-fu ist nicht bis zur Aufgabe des Übersetzens dieses Wissens in der effizienten Code
  2. Während ich eine binomische Kurve würde ein genaues Ergebnis zu verstehen, habe ich den Eindruck, dass sie von Natur aus ineffizient wäre. Eine schnelle Methode, um ein ungefähres Ergebnis berechnen würde genügen.

Ich möchte betonen, dass diese Berechnung muss schnell sein und sollte idealerweise mit Standard 64 oder 128-Bit-Gleitkomma-Berechnung bestimmbar sein.

Ich suche nach einer Funktion, die P, S und N nimmt - und gibt eine Wahrscheinlichkeit. Da ich mehr vertraut mit dem Code bin als mathematische Notation, würde ich es vorziehen, dass alle Antworten Pseudo-Code oder Code verwenden.

War es hilfreich?

Lösung

Exact Binomialverteilung 

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

Normale Schätzung, 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: Gut 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

Andere Tipps

Ich war an einem Projekt, wo wir in der Lage sein notwendig, um die binomische CDF in einer Umgebung zu berechnen, die keine faktoriellen oder Gamma-Funktion definiert hatte. Es dauerte ein paar Wochen, aber ich landete mit dem folgenden Algorithmus kommen, die die CDF genau berechnet (das heißt keine Angleichung erforderlich). Python ist im Grunde so gut wie Pseudo-Code, nicht wahr? 

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

Performance 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 will nicht in eine vollständige Ableitung dieses Algorithmus gehen, weil nicht Stackoverflow Mathjax nicht unterstützt, aber die Stoßrichtung wird zunächst die folgende Äquivalenz zu identifizieren:

  • Für alle k> 0, sp.misc.comb(n,k) == np.prod([(n-k+1)/k for k in range(1,k+1)])

Was können wir umschreiben als:

  • sp.misc.comb(n,k) == sp.misc.comb(n,k-1) * (n-k+1)/k

oder in Protokollspeicher:

  • np.log( sp.misc.comb(n,k) ) == np.log(sp.misc.comb(n,k-1)) + np.log(n-k+1) - np.log(k)

Da die CDF eine Summierung PMFs ist, können wir diese Formulierung verwenden, um die binomischen Koeffizienten (das Protokoll von denen b in der obigen Funktion ist) zu berechnen für PMF_ {x = i} von dem Koeffizienten wir für PMF_ berechnet {x = i-1}. Dies bedeutet, dass wir alles in einer einzigen Schleife mit Akkumulatoren tun können, und wir brauchen keine factorials zu berechnen!

Der Grund, die meisten der Berechnungen in Log-Raum durchgeführt werden, ist die numerische Stabilität der Polynomtermen zu verbessern, das heißt p^x und (1-p)^(1-x) das Potenzial haben, sehr groß oder sehr klein zu sein, was Rechenfehler verursachen kann.

EDIT: Ist dies ein neuer Algorithmus? Ich habe seit etwa ein und aus Stossen, bevor ich gepostet, und ich frage mich immer, wenn ich dies auf mehr formell schreiben soll und lege ihn auf eine Zeitschrift.

Ich glaube, Sie die unvollständige Beta-Funktion auswerten möchten.

Es gibt eine schöne Umsetzung eine Kettenbruch Darstellung in "Numerical Recipes in C" verwenden, Kapitel 6:. 'Sonderfunktionen'

Ich kann nicht völlig bürgt für die Effizienz, sondern Scipy hat ein Modul für diese 

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

Ein effiziente und, was noch wichtiger ist, numerisch stabiler Algorithmus existiert im Bereich der Bezier-Kurven verwendet in Computer Aided Design. Es heißt de Casteljau-Algorithmus verwendet die Bernstein Polynomials zur Bewertung Bézier-Kurven definieren.

Ich glaube, dass ich nur pro Antwort einen Link erlaubt bin so beginnen mit Wikipedia - Bernstein Polynomials

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

  

Lets sagen, ich weiß, dass die Wahrscheinlichkeit, ein Heads mit einer bestimmten Münze zu werfen, ist P.   Wie hoch ist die Wahrscheinlichkeit, mich zu werfen   die Münze T-mal und immer mindestens   S Köpfe?

  • Set n = T
  • Set Beta [i] = 0 für i = 0, ... S - 1
  • Set Beta [i] = 1 für i = S, ... T
  • Stellen t = p
  • Ausrechnen B (t) unter Verwendung von de Casteljau
  

oder höchstens S Köpfe?

  • Set n = T
  • Set Beta [i] = 1 für i = 0, ... S
  • Set Beta [i] = 0 für i = S + 1, ... T
  • Stellen t = p
  • Ausrechnen B (t) unter Verwendung von de Casteljau

Open-Source-Code wahrscheinlich ist bereits vorhanden. NURBS-Kurven (Non-Uniform Rational B-Spline-Kurven) ist eine Verallgemeinerung von Bezier-Kurven und wird in CAD weit verbreitet. Versuchen Sie openNURBS (die Lizenz ist sehr liberal) oder andernfalls Open CASCADE (etwas weniger liberal und undurchsichtig Lizenz). Beide Toolkits sind in C ++, aber IIRC, .NET Bindungen existieren.

Wenn Sie Python verwenden, keine Notwendigkeit, es selbst zu codieren. Scipy hat für Sie: 

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

Aus dem Teil Ihrer Frage „immer mindestens S Köpfe“ wollen Sie die kumulative Binomialverteilung Funktion. Siehe http://en.wikipedia.org/wiki/Binomial_distribution für die Gleichung, die ist, beschrieben, wie in Bezug auf die „regularisiert unvollständigen beta-Funktion“ (wie bereits beantwortete) ist. Wenn Sie nur die Antwort berechnen wollen, ohne die gesamte Lösung selbst, die GNU Scientific Library bietet die Funktion zu implementieren. Gsl_cdf_binomial_P und gsl_cdf_binomial_Q

Das DCDFLIB Projekt C # Funktionen (Wrapper um C-Code) hat viele zu bewerten CDF-Funktionen, einschließlich der Binomialverteilung. Sie können das Original C und Fortran-Code finden hier . Dieser Code ist gut getestet und genau sind.

Wenn Sie Ihren eigenen Code schreiben möchten, auf eine externe Bibliothek abhängig ist zu vermeiden, könnte man die normale Annäherung an die in anderen Antworten erwähnt binomischen verwenden. Hier sind einige Hinweise auf , wie gut die Annäherung unter verschiedenen Umständen wird. Wenn Sie diesen Weg gehen und Code, um die normalen CDF zu berechnen, ist hier Python-Code dafür, dass . Es ist nur etwa ein Dutzend Zeilen Code und leicht an jede andere Sprache portiert werden könnten. Aber wenn Sie eine hohe Genauigkeit und effizienten Code wollen, sind Sie besser dran, Code von Drittanbietern wie DCDFLIB verwenden. Mehrere Mannjahre ging in die Bibliothek zu erzeugen.

Versuchen Sie diese , verwendet in GMP. Eine weitere Referenz ist diese . 

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.
Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top