Como posso de forma eficiente calcular a função de distribuição cumulativa binomial?
-
11-09-2019 - |
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:
- A minha matemática-fu não está à altura da tarefa de traduzir esse conhecimento em código eficiente
- 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.
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.
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.