Como implementar a função hipergeométrica generalizada para usar em beta-binomial cdf, sf, ppf?

7

Estou escrevendo uma subclasse scipy.stats._distn_infrastructure.rv_discretepara a distribuição binomial beta cujo PMF é

P(X=kN,α,β)(Nk)B(k+α,Nk+β)B(α,β),

onde B é a função Beta. Minha implementação atual do CDF e SF (função de sobrevivência, equivalente a 1 - CDF) é imprecisa; a estratégia que empreguei calcula o valor esperado do binomial cdf em relação ao componente beta:

PBB(XkN,α,β)=Ep[PBinom(XkN,p)],
que pBeta(α,β) . Consigo isso usando o scipy.stats.beta.expectmétodo, que não é vetorizado de forma inata (ele trava em qualquer coisa que não seja uma matriz float ou 0d).

O PPF é ainda pior - é um loop de força bruta sobre os números inteiros k=0,,N tal que

P(XkN,α,β)q.

Segundo a Wikipedia, a função de sobrevivência para a distribuição beta-binomial é

P(X>kN,α,β)=B(β+nk1,α+k+1)3F2(a,b;k)B(α,β)B(nk,k+2)(n+1),

onde é a função hipergeométrica generalizada. Existe uma maneira eficiente de calcular isso em Python, para que eu possa remover a referência ? Além disso, como eu inverteria essa função para resolver dado ?3F2beta.expectkq=P(XkN,α,β)

Scott Norton
fonte
Talvez seja útil saber que, para os valores de que (implicitamente) aparecem aqui, é um polinômio em (de grau , ). Não simplifica em geral. a,b3F2(;;z)znk11kn1
whuber
Você encontrou alguma solução para sua pergunta? Se sim, talvez você queira compartilhá-lo como resposta à sua pergunta?
Tim

Respostas:

2

Isso não responde diretamente à sua pergunta, mas se você estiver pensando em estimar a função de distribuição cumulativa do binômio beta com mais eficiência, poderá usar um algoritmo recursivo que é um pouco mais eficiente que a implementação ingênua.

Observe que a função massa de probabilidade da distribuição beta-binomial

f(x)=(nx)B(x+α,nx+β)B(α,β)

pode ser reescrito se você lembrar que ee que , para que se torneB(x,y)=Γ(x)Γ(y)Γ(x+y)Γ(x)=(x1)!(nk)=i=1kn+1ii

f(x)=(i=1xn+1ii)(α+x1)!(β+nx1)!(α+β+n1)!B(α,β)

Isso faz com que a atualização de para fácilxx+1

f(x+1)=(i=1xn+1ii)n+1x+1x+1(α+x1)!(α+x)(β+nx1)!(β+nx)1(α+β+n1)!(α+β+n)B(α,β)

e usando isso, você pode calcular a função de distribuição cumulativa como

F(x)=k=0xf(k)

usando apenas operações aritméticas simples, em vez de calcular funções mais intensivas em computador.

Sidenote: ao lidar com grandes números, você entraria em problemas de precisão numérica, portanto, um código mais robusto precisaria trabalhar com logaritmos, mas mesmo que você esperasse uma melhoria na eficiência (código até duas a três vezes mais rápido quando eu executei alguns benchmarks em Código C ++ implementando-o em comparação com a implementação ingênua).

Tim
fonte
11
Outra nota. A razão de integrais beta para o primeiro termo é outro produto simples isso simplifica paraf(0)=B(a,n+b)B(a,b)=Γ(n+b)Γ(a+b)Γ(n+a+b)Γ(b)j=1nn+bjn+a+bj
probabilityislogic