Come posso calcolare in modo efficiente la funzione di distribuzione cumulativa binomiale?

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

Domanda

Diciamo che so che la probabilità di un "successo" è P.Eseguo il test N volte e vedo S successi.Il test è simile al lancio di una moneta dal peso non uniforme (forse testa è un successo, croce è un fallimento).

Voglio conoscere la probabilità approssimativa di vedere S successi o un numero di successi meno probabili di S successi.

Quindi, ad esempio, se P è 0,3, N è 100 e ottengo 20 successi, cerco la probabilità di ottenerne 20 o meno successi.

Se, d'altro canto, P è 0,3, N è 100 e ottengo 40 successi, cerco la probabilità di ottenere 40 nostri successi in più.

Sono consapevole che questo problema riguarda la ricerca dell'area sotto una curva binomiale, tuttavia:

  1. Il mio math-fu non è all’altezza del compito di tradurre questa conoscenza in un codice efficiente
  2. Anche se capisco che una curva binomiale darebbe un risultato esatto, ho l'impressione che sarebbe intrinsecamente inefficiente.Sarebbe sufficiente un metodo veloce per calcolare un risultato approssimativo.

Dovrei sottolineare che questo calcolo deve essere veloce e idealmente dovrebbe essere determinabile con il calcolo in virgola mobile standard a 64 o 128 bit.

Sto cercando una funzione che prenda P, S e N e restituisca una probabilità.Dato che ho più familiarità con il codice che con la notazione matematica, preferirei che qualsiasi risposta utilizzi pseudo-codice o codice.

È stato utile?

Soluzione

esatta distribuzione binomiale 

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 Stima, buono per n grande 

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 Stima: Per n grande e piccola 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

Altri suggerimenti

ero su un progetto in cui abbiamo bisogno di essere in grado di calcolare il binomio CDF in un ambiente che non ha avuto una funzione fattoriale o gamma definita. Mi ci sono voluti un paio di settimane, ma ho finito per venire con il seguente algoritmo che calcola la CDF esattamente (vale a dire senza approssimazione necessario). Python è fondamentalmente buono come pseudocodice, giusto? 

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 scale con x. Per piccoli valori di x, questa soluzione è di circa un ordine di grandezza più veloce di scipy.stats.binom.cdf, con prestazioni simili a circa x = 10.000.

Non voglio entrare in una derivazione completa di questo algoritmo perché StackOverflow non supporta mathjax, ma la spinta di esso è il primo che identifica la seguente equivalenza:

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

Il che possiamo riscrivere come:

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

o in spazio di log:

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

Poiché il CDF è un riepilogo di PMF, possiamo usare questa formulazione per calcolare il coefficiente binomiale (il registro dei quali è b nella funzione sopra) per PMF_ {x = i} dal coefficiente abbiamo calcolato per PMF_ {x = i-1}. Questo significa che possiamo fare tutto all'interno di un unico ciclo con accumulatori, e non abbiamo bisogno di calcolare eventuali fattoriali!

Il motivo più dei calcoli sono fatti in spazio di registrazione è di migliorare la stabilità numerica dei termini polinomiali, cioè p^x e (1-p)^(1-x) hanno il potenziale per essere estremamente grande o estremamente piccolo, che può causare errori di calcolo.

EDIT: Si tratta di un nuovo algoritmo? Sto rovistando dentro e fuori da prima ho postato questo, e sto sempre chiedo se devo scrivere questo in modo più formale e la sottopone ad un giornale.

Penso che si desidera valutare le href="http://mathworld.wolfram.com/BinomialDistribution.html" incompleta della funzione beta .

C'è una bella applicazione utilizzando una rappresentazione frazione continua in "Numerical Recipes In C", capitolo 6:. 'Funzioni speciali'

Non posso assolutamente garantire per l'efficienza, ma SciPy ha un modulo per questo 

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

Esiste un algoritmo efficiente e, soprattutto, numericamente stabile nel dominio di Curve di Bezier utilizzato nella progettazione assistita da computer.È chiamato Algoritmo di de Casteljau utilizzato per valutare il Polinomi di Bernstein utilizzato per definire le curve di Bezier.

Credo che mi sia consentito un solo collegamento per risposta, quindi inizia con Wikipedia - Polinomi di Bernstein

Notare la relazione molto stretta tra la distribuzione binomiale e i polinomi di Bernstein.Quindi fare clic sul collegamento sull'algoritmo di de Casteljau.

Diciamo che so che la probabilità di ottenere testa con una particolare moneta è P.Qual è la probabilità che io lanci i tempi di moneta e mi ottenga almeno la testa?

  • Imposta n = T
  • Imposta beta[i] = 0 per i = 0, ...S-1
  • Imposta beta[i] = 1 per i = S, ...T
  • Imposta t = p
  • Valutare B(t) utilizzando de Casteljau

o al massimo S teste?

  • Imposta n = T
  • Imposta beta[i] = 1 per i = 0, ...S
  • Poniamo beta[i] = 0 per i = S + 1, ...T
  • Imposta t = p
  • Valutare B(t) utilizzando de Casteljau

Probabilmente il codice open source esiste già. Curve NURBS (Curve B-spline razionali non uniformi) sono una generalizzazione delle curve di Bezier e sono ampiamente utilizzate nel CAD.Prova openNurbs (la licenza è molto liberale) o, in mancanza, Open CASCADE (una licenza un po' meno liberale e opaca).Entrambi i toolkit sono in C++, tuttavia esistono collegamenti IIRC e .NET.

Se si sta utilizzando Python, nessuna necessità di codificare da soli. SciPy abbiamo coperto: 

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

Dalla parte della tua domanda "ottenere almeno teste S" si desidera che la funzione di distribuzione binomiale cummulative. Vedere http://en.wikipedia.org/wiki/Binomial_distribution per l'equazione, che è descritto come essere in termini di "regolarizzato funzione beta incompleta" (come già risposto). Se si desidera solo per calcolare la risposta senza dover implementare l'intera soluzione da soli, la Biblioteca Scientifica GNU fornisce la funzione:. Gsl_cdf_binomial_P e gsl_cdf_binomial_Q

Il DCDFLIB Progetto ha funzioni # C (involucri intorno codice C) per valutare molti funzioni CDF, tra cui la distribuzione binomiale. È possibile trovare la C originale e il codice FORTRAN qui . Questo codice è ben collaudato e preciso.

Se si desidera scrivere il proprio codice per evitare di essere dipendente da una libreria esterna, è possibile utilizzare l'approssimazione normale al binomio detto in altre risposte. Ecco alcune note su quanto è buono il ravvicinamento è in varie circostanze. Se seguire questa strada e hanno bisogno di codice per calcolare il normale CDF, ecco Python codice per farlo . E 'solo una dozzina di righe di codice e potrebbe essere facilmente portato su qualsiasi altra lingua. Ma se si vuole elevata precisione e codice efficiente, è meglio utilizzare il codice di terze parti come DCDFLIB. Diversi anni-uomo è andato in produzione di quella libreria.

questo , utilizzato in GMP. Un altro riferimento è questo . 

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.
Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top