Como posso de forma eficiente calcular a função de distribuição cumulativa binomial?

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

Pergunta

Vamos dizer que eu sei que a probabilidade de um "sucesso" é P. eu corro os tempos de teste N, e vejo sucessos S. O teste é semelhante a jogar uma moeda desigualmente ponderada (talvez cabeças é um sucesso, caudas é um fracasso).

Eu quero saber a probabilidade aproximada de ver tanto sucessos S, ou um número de sucessos menos provável do que sucessos S.

Assim, por exemplo, se P é de 0,3, N é de 100, e fico com 20 sucessos, eu estou procurando a probabilidade de obter 20 ou menos sucessos.

Se, por outro hadn, P é de 0,3, N é de 100, e fico com 40 sucessos, eu estou procurando a probabilidade de obter 40 nossos mais sucessos.

Estou ciente de que este problema diz respeito a encontrar a área sob uma curva binomial, no entanto:

  1. A minha matemática-fu não está à altura da tarefa de traduzir esse conhecimento em código eficiente
  2. Embora compreenda uma curva binomial daria um resultado exato, tenho a impressão de que seria inerentemente ineficiente. Um método rápido para calcular um resultado aproximado seria suficiente.

I deve salientar que este cálculo tem de ser rápido, e idealmente deve ser determinável com a computação ponto flutuante padrão de 64 ou 128 bits.

Eu estou procurando uma função que leva P, S e N - e retorna uma probabilidade. Como eu estou mais familiarizado com código de notação matemática, eu preferiria que qualquer pseudo-código respostas empregar ou código.

Foi útil?

Solução

Exact Distribuição 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 Estimativa, bom para grande 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 Estimado: Bom para n grande e pequena 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

Outras dicas

Eu estava em um projeto onde nós precisava de ser capaz de calcular o CDF binomial em um ambiente que não tem uma função fatorial ou gama definida. Levei algumas semanas, mas acabei chegando com o seguinte algoritmo que calcula a CDF exatamente (ou seja, não necessária aproximação). Python é basicamente tão bom quanto pseudocódigo, certo? 

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

escalas de desempenho com x. Para valores pequenos de x, esta solução é de cerca de uma ordem de magnitude mais rápido do que scipy.stats.binom.cdf, com desempenho semelhante em torno de x = 10.000.

Eu não vou entrar em uma derivação completa deste algoritmo porque stackoverflow não suporta MathJax, mas o impulso do que é primeiro identificar o seguinte equivalência:

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

O que podemos reescrever como:

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

ou no espaço de log:

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

Como a CDF é um somatório de PMF, podemos utilizar esta formulação para calcular o coeficiente binomial (o registo dos quais é b na função acima) para PMF_ {x = i} a partir do coeficiente de nós calculada para PMF_ {x = i-1}. Isto significa que podemos fazer tudo dentro de um único loop usando acumuladores, e nós não precisa calcular quaisquer fatoriais!

A razão a maioria dos cálculos são feitos no espaço de log é melhorar a estabilidade numérica dos termos polinomiais, ou seja p^x e (1-p)^(1-x) tem o potencial para ser muito grande ou muito pequeno, o que pode causar erros de cálculo.

EDIT: Este é um novo algoritmo? Eu tenho picar ao redor dentro e fora desde antes de eu postei isso, e eu estou cada vez mais querendo saber se eu deveria escrever isso mais formalmente e submetê-lo a uma revista.

Eu acho que você quer avaliar a incompleta função beta .

Há uma implementação agradável usando uma representação fração contínua em "Numerical Recipes In C", capítulo 6:. 'Funções especiais'

Não consigo totalmente atestar a eficiência, mas Scipy tem um módulo para este 

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

Um eficiente e, mais importante, o algoritmo estável numérica existe no domínio da Curvas Bezier usado em Desenho Assistido por Computador. Ele é chamado de algoritmo de Casteljau usado para avaliar a Bernstein Polinômios usado para definir curvas de Bezier.

Acredito que estou apenas permitido um link por resposta assim que começar com Wikipedia - Bernstein Polinômios

Observe a relação muito estreita entre a distribuição binomial e polinômios Bernstein. Em seguida, clicar no link em de algoritmo de Casteljau.

Vamos dizer que eu sei que a probabilidade de jogar um heads com uma moeda em particular é P. Qual é a probabilidade de me atirar os tempos moeda T e recebendo pelo menos S cabeças?

  • Set n = T
  • Set beta [i] = 0 para i = 0, ... S - 1
  • Set beta [i] = 1 para i = S, ... T
  • Set t = p
  • Avaliar B (t) usando de Casteljau

ou na maioria das cabeças S?

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

Abrir código-fonte, provavelmente já existe. NURBS Curves (Rational B-spline Curves não uniforme) são uma generalização das curvas de Bezier e são amplamente utilizados em CAD. Tente OpenNURBS (a licença é muito liberal) ou deixar que o Open CASCADE (a licença um pouco menos liberal e opaco). Ambos os conjuntos de ferramentas são em C ++, embora, IIRC, NET existem ligações.

Se você estiver usando Python, não há necessidade de código-lo sozinho. Scipy temos que cobertas: 

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 da parte da sua pergunta "chegando cabeças S menos" você quer a função de distribuição binomial cumulativa. Consulte http://en.wikipedia.org/wiki/Binomial_distribution para a equação, que é descrita como sendo em termos da "regularizado função beta incompleta" (como já respondeu). Se você quiser apenas para calcular a resposta sem ter de implementar toda a solução sozinho, o GNU Scientific Library fornece a função:. Gsl_cdf_binomial_P e gsl_cdf_binomial_Q

O DCDFLIB Projeto tem funções C # (invólucros em torno do código C) para avaliar muitos CDF funções, incluindo a distribuição binomial. Você pode encontrar o C original e código Fortran aqui . Este código é bem testado e precisa.

Se você quiser escrever seu próprio código para evitar ser dependente de uma biblioteca externa, você pode usar a aproximação normal para a binomial mencionado em outras respostas. Aqui estão algumas notas sobre quão boa a aproximação é sob várias circunstâncias. Se você ir por esse código de rota e necessidade de calcular o CDF normal, aqui está Python código para fazer isso . É apenas cerca de uma dúzia de linhas de código e pode ser facilmente transportado para qualquer outra língua. Mas se você quiser alta precisão e código eficiente, você é melhor fora de usar código de terceiros, como DCDFLIB. Vários anos-homem entrou em produzir essa biblioteca.

Tente esta , usado em GMP. Outra referência é 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 em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top