¿Cómo puedo calcular de manera eficiente la función de distribución acumulativa binomial?

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

Pregunta

Digamos que sé la probabilidad de un "éxito" es P. funciono con la prueba de N veces, y veo S éxitos. La prueba es similar a lanzar una moneda de manera desigual ponderada (quizás cabezas es un éxito, colas es un fracaso).

Quiero saber la probabilidad aproximada de ver cualquiera de éxitos S, o una serie de éxitos menos probabilidades que los éxitos S.

Así, por ejemplo, si P es 0,3, N es 100, y me da 20 éxitos, estoy buscando la probabilidad de obtener 20 o menos éxitos.

Si, por otro hadn, P es 0,3, N es 100, y me da 40 éxitos, estoy buscando la probabilidad de obtener 40 nuestros más éxitos.

Soy consciente de que este problema se refiere a encontrar el área bajo una curva binomial, sin embargo:

  1. Mi matemáticas-fu no es hasta la tarea de traducir este conocimiento en código eficiente
  2. Si bien entiendo una curva binomial daría un resultado exacto, me da la impresión de que sería inherentemente ineficiente. Un método rápido para calcular un resultado aproximado sería suficiente.

he de destacar que este cálculo tiene que ser rápido, e idealmente debería ser determinable con el cálculo estándar punto 64 o 128 bits flotante.

Busco una función que toma P, S, y N - y devuelve una probabilidad. Como estoy más familiarizado con el código de notación matemática, preferiría que ninguna respuesta emplean pseudo-código o código.

¿Fue útil?

Solución

Exact Distribución Binomial 

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

Normal Estimación, buena para 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 Estimado: Bueno para n grande y pequeña 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

Otros consejos

Yo estaba en un proyecto en el que teníamos que estar en condiciones de calcular el binomio CDF en un ambiente que no tienen una función factorial o gamma definido. Me tomó un par de semanas, pero terminó por dar con el siguiente algoritmo que calcula la CDF exactamente (es decir, no necesaria aproximación). Python es básicamente tan bueno como pseudocódigo, ¿verdad? 

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

Rendimiento escala con x. Para valores pequeños de x, esta solución es aproximadamente un orden de magnitud más rápido que scipy.stats.binom.cdf, con un rendimiento similar en torno a x = 10.000.

No voy a entrar en una derivación completa de este algoritmo porque stackoverflow no admite mathjax, pero la idea central de que es identificar primero la siguiente equivalencia:

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

¿Qué se puede reescribir como:

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

o en el espacio de registro:

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

Debido a que el CDF es una suma de FPC, podemos utilizar esta formulación para calcular el coeficiente binomial (el registro de las cuales es b en la función anterior) para PMF_ {x = i} a partir del coeficiente se calculó para PMF_ {x = i-1}. Esto significa que podemos hacer todo dentro de un solo bucle mediante acumuladores, y no necesitamos para calcular los factoriales!

La razón por la mayoría de los cálculos se realizan en el espacio de registro es mejorar la estabilidad numérica de los términos del polinomio, es decir p^x y (1-p)^(1-x) tienen el potencial de ser muy grande o muy pequeña, lo que puede causar errores de cálculo.

EDIT: ¿Es este un nuevo algoritmo? He estado hurgando de forma intermitente desde antes he publicado esto, y me pregunto si cada vez que debería escribir esto de manera más formal y enviarlo a una revista.

Creo que desea evaluar los href="http://mathworld.wolfram.com/BinomialDistribution.html" función beta incompleta .

Hay una buena aplicación utilizando una representación en fracción continua "Numerical Recipes En C", capítulo 6:. 'Funciones especiales'

No puedo responder por completo por la eficiencia, pero Scipy tiene un módulo para este 

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

Una eficiente y, más importante aún, existe algoritmo estable numérica en el campo de Curvas de Bezier utilizan en diseño asistido por ordenador. Se llama de algoritmo de Casteljau Se utiliza para evaluar la Bernstein polinomios se utiliza para definir las curvas de Bézier.

Creo que sólo se me permite un enlace por respuesta para comenzar con Wikipedia - Bernstein polinomios

Note la relación muy estrecha entre la distribución binomial y la Bernstein polinomios. A continuación, hacer clic a través del enlace en el algoritmo de de Casteljau.

  

Digamos que sé la probabilidad de que salga un mano a mano con una moneda en particular es P.   ¿Cuál es la probabilidad de lanzándome   los tiempos de la moneda T y conseguir, al menos,   S cabezas?

  • Set n = T
  • beta Ajuste [i] = 0 para i = 0, ... S - 1
  • beta Ajuste [i] = 1 para i = S, ... T
  • Establecer T = p
  • Evaluar B (t) utilizando de Casteljau
  

o en la mayoría de los jefes S?

  • Set n = T
  • beta Ajuste [i] = 1 para i = 0, ... S
  • Set beta [i] = 0 para i = S + 1, ... T
  • Establecer T = p
  • Evaluar B (t) utilizando de Casteljau

código fuente abierto, probablemente ya existe. curvas NURBS (racionales no uniformes curvas B-spline) son una generalización de las curvas de Bézier y son ampliamente utilizados en CAD. Trate openNURBS (la licencia es muy liberal) o en su defecto abierto en cascada (una licencia algo menos liberal y opaco). Ambas herramientas están en C ++, sin embargo, existen fijaciones IIRC, .NET.

Si está utilizando Python, sin necesidad de código usted mismo. Scipy tiene todo cubierto: 

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

A partir de la parte de su pregunta "conseguir al menos cabezas S" desea que la función de distribución acumulativa binomial. Ver http://en.wikipedia.org/wiki/Binomial_distribution para la ecuación, la cual es descrito como en términos de la "función beta incompleta regularizada" (como ya se ha contestado). Si lo que desea es calcular la respuesta sin tener que implementar toda la solución a sí mismo, la Biblioteca GNU Scientific ofrece la función:. Gsl_cdf_binomial_P y gsl_cdf_binomial_Q

El DCDFLIB Proyecto tiene funciones de C # (envolturas alrededor de código C) para evaluar muchos funciones CDF, incluyendo la distribución binomial. Usted puede encontrar el original de C y FORTRAN aquí . Este código está bien probado y precisa.

Si desea escribir su propio código para evitar depender de una biblioteca externa, se puede usar la aproximación normal a la binomial se ha mencionado en otras respuestas. Aquí hay algunas notas sobre lo bien que la aproximación se en diversas circunstancias. Si ir por ese camino y necesita el código para calcular la FDA normal, aquí está Python código para hacer eso . Es sólo una docena de líneas de código y podría fácilmente ser portado a cualquier otro idioma. Pero si quieres una alta precisión y un código eficiente, es mejor utilizar código de terceros como DCDFLIB. Varios años-hombre entró en la producción de esa biblioteca.

Trate éste , utilizado en GMP. Otra referencia es este . 

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.
Licenciado bajo: CC-BY-SA con atribución
No afiliado a StackOverflow
scroll top