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

Question

Digamos que sé que la probabilidad de un "éxito" es P. Ejecuto la prueba N veces, y veo S éxitos. La prueba es similar a arrojar una moneda desigualmente ponderada (quizás las cabezas son un éxito, las colas son un fracaso).

Quiero saber la probabilidad aproximada de ver los éxitos S, o un número de éxitos menos probable que los éxitos S. Por ejemplo, si P es 0.3, N es 100 y obtengo 20 éxitos, estoy buscando la probabilidad de obtener 20 o menos éxitos en.

Si, por otro lado, P es 0,3, N es 100 y obtengo 40 éxitos, estoy buscando la probabilidad de obtener 40 más éxitos. Sin embargo

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

1. Mi matemáticas-fu no es hasta la tarea de traducir este conocimiento en código eficiente
2. Mientras Entiendo que 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.

Debo recalcar que este cálculo tiene que ser rápido, y idealmente podría determinarse con un cálculo de coma flotante de 64 o 128 bits estándar.

Estoy buscando una función que tome P, S y N, y devuelve una probabilidad. Como estoy más familiarizado con el código que con la notación matemática, preferiría que cualquier respuesta emplee pseudo-código o código.

Answer 1

exacta 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 estimado, 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 grandes valores de n y p pequeño

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 

Answer 2

Probar this one, utilizado en GMP. Otra referencia es this.

Answer 3

Creo que quiere evaluar el incomplete beta function.

Hay una buena implementación usando una representación de fracción continua en "Recetas numéricas en C", capítulo 6: "Funciones especiales".

Answer 4

De la parte de su pregunta "obtener al menos S cabezas" desea la función de distribución binomial cumulativa. Ver http://en.wikipedia.org/wiki/Binomial_distribution para la ecuación, que se describe como en términos de la "función beta incompleta regularizada" (como ya se respondió). Si solo desea calcular la respuesta sin tener que implementar la solución completa usted mismo, la Biblioteca Científica GNU proporciona la función: gsl_cdf_binomial_P y gsl_cdf_binomial_Q.

Answer 5

Existe un algoritmo numérico eficiente y, más importante aún, numérico en el dominio de Bezier Curves utilizado en Computer Aided Design. Se llama del algoritmo de Casteljau utilizado para evaluar Bernstein Polynomials utilizado para definir Curvas de Bezier.

creo que sólo se me permite un enlace por respuesta así que empieza con Wikipedia - Bernstein Polynomials

Aviso la estrecha relación entre la distribución binomial y el Bernstein polinomios. Luego haz clic en el enlace del algoritmo de Casteljau.

Digamos que sé que 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 T moneda y obtener al menos cabezas S?

• conjunto N = T
• Conjunto beta [i] = 0 para i = 0, ... S - 1
• Conjunto beta [i] = 1 para i = S, .. . T
• conjunto T = p
• Evaluar B (t) utilizando de Casteljau

o en la mayoría de los jefes S?

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

código fuente abierto, probablemente ya existe. Curvas NURBS (Curvas B-spline racionales no uniformes) son una generalización de Curvas de Bezier y son ampliamente utilizadas en CAD. Pruebe openNurbs (la licencia es muy liberal) o, en su defecto, Open CASCADE (una licencia algo menos liberal y opaca). Ambos kits de herramientas están en C++, sin embargo, existen enlaces IIRC, .NET.

Answer 6

El DCDFLIB Project tiene funciones C# (envoltorios alrededor del código C) para evaluar muchas funciones de CDF, incluida la distribución binomial. Puede encontrar el código original de C y FORTRAN here. Este código está bien probado y es preciso.

Si desea escribir su propio código para evitar depender de una biblioteca externa, puede usar la aproximación normal al binomio mencionado en otras respuestas. Aquí hay algunas notas en how good the approximation is bajo diversas circunstancias. Si sigue esa ruta y necesita un código para calcular el CDF normal, aquí está Python code por hacerlo. Solo se trata de una docena de líneas de código y podría portarse fácilmente a cualquier otro idioma. Pero si desea una alta precisión y un código eficiente, es mejor que utilice un código de terceros como DCDFLIB. Varios años-hombre entraron en la producción de esa biblioteca.

Answer 7

No puedo responder por completo por la eficiencia, pero tiene un Scipy module for this

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

Answer 8

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 

Answer 9

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ó algunas semanas, pero terminé con el siguiente algoritmo que calcula el CDF exactamente (es decir, no es necesaria una aproximación). Python es básicamente tan bueno como el 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 

Balanzas de rendimiento con x. Para valores pequeños de x, esta solución tiene 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 soporta mathjax, pero la idea central de que es identificar primero la siguiente equivalencia:

• para todo 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 usar esta formulación para calcular el coeficiente binomial (el registro de las cuales es b en el función arriba) para PMF_ {x = i} del coeficiente que calculamos para PMF_ {x = i-1}. Esto significa que podemos hacer todo en un solo bucle utilizando acumuladores, ¡y no necesitamos calcular ningún factorial! El motivo por el que la mayoría de los cálculos se realizan en el espacio de registro es mejorar la estabilidad numérica de los términos polinomiales, es decir, p^x y (1-p)^(1-x) tienen el potencial de ser extremadamente grandes o extremadamente pequeños, lo que puede provocar errores de cálculo.