Como amostrar da distribuição discreta nos números inteiros não negativos?

9

Eu tenho a seguinte distribuição discreta, onde são constantes conhecidas:α,β

p(x;α,β)=Beta(α+1,β+x)Beta(α,β)for x=0,1,2,

Quais são algumas das abordagens para obter amostras eficientes dessa distribuição?

jII
fonte

Respostas:

9

Esta é uma distribuição binomial beta negativa , com o parâmetro no seu caso, usando a notação da Wikipedia. Também denominou distribuição Beta-Pascal quando é um número inteiro. Como você observou em um comentário, esta é uma distribuição preditiva no modelo binomial negativo bayesiano com um Beta conjugado anterior à probabilidade de sucesso.r=1r

Assim, você pode amostrá-lo amostrando uma variável e, em seguida, amostrando uma variável binomial negativa (com no seu caso, ou seja, para dizer uma distribuição geométrica).Beta(α,β)uNB(r,u)r=1

Esta distribuição é implementada no pacote de R brr. O amostrador tem nome rbeta_nbinom, pmf tem nome dbeta_nbinom, etc. As notações são , , . Verifica:a=rc=αd=β

> Alpha <- 2; Beta <- 3
> a <- 1
> all.equal(brr::dbeta_nbinom(0:10, a, Alpha, Beta), beta(Alpha+a, Beta+0:10)/beta(Alpha,Beta))
[1] TRUE

Observando o código, podemos ver que ele realmente chama a ghyperfamília de distribuições (hipergeométricas generalizadas) do SuppDistspacote:

brr::rbeta_nbinom
function(n, a, c, d){
  rghyper(n, -d, -a, c-1)
}

De fato, a distribuição do BNB é conhecida como distribuição hipergeométrica generalizada do tipo IV . Veja a ajuda de ghyperno SuppDistspacote. Acredito que isso também pode ser encontrado no livro da Johnson & Al, Univariate Discrete Distributions .

Stéphane Laurent
fonte
Essa resposta é ótima, mas seria ainda melhor se você provar que o OP de densidade postado é o mesmo que a densidade binomial negativa.
Sycorax diz Restabelecer Monica
11
@ user777 Acho que o autor do OP provou a si próprio, tendo em vista seu comentário à resposta de Xian (distribuição preditiva posterior no modelo binomial negativo com um beta anterior conjugado).
Stéphane Laurent
10

Dado que está diminuindo com , sugiro gerar uma variável uniforme e calculando as somas acumuladas até A realização é então igual ao correspondente . Desde a

Beta(α+1,β+x)Beta(α,β)=αα+β+xβ+x1α+β+x1βα+β
xuU(0,1)
Sk=x=0kBeta(α+1,β+x)Beta(α,β)
Sk>u
k
Rx=Beta(α+1,β+x)Beta(α,β)=αα+β+xβ+x1α+β+x1βα+β=α+β+x1α+β+xβ+x1α+β+x1Rx1=β+x1α+β+xRx1
e o cálculo pode evitar o uso completo das funções Gamma.
Sk=Sk1+Rk
Xi'an
fonte
11
(+1) Usar agilizará bastante o trabalho. Sk=1Γ(a+b)Γ(b+k+1)Γ(b)Γ(a+b+k+1)
whuber
11
Re edição: Eu suspeito que explorar as funções Gamma, no entanto, será útil para resolver em termos de , e . Por exemplo, pode-se encontrar uma aproximação inicial de usando a fórmula de Stirling na avaliação de e e depois polindo-a com alguns Newton-Raphson passos. Aqueles precisam avaliar o log Gamma e seus derivados. É claro que se e são integrais, então a solução é a raiz de um polinômio - mas mesmo assim, usar Gamma ainda pode ser o caminho a percorrer. u α β u Γ ( b + k + 1 ) Γ ( a + b + k + 1 ) α βkuαβuΓ(b+k+1)Γ(a+b+k+1)αβ
whuber
11
Ótima resposta! Aceitei a resposta dada pelo SL porque trouxe à minha atenção um ponto-chave (não faz parte da pergunta original), que a amostragem a partir do preditivo posterior é equivalente a amostrar o parâmetro a partir do posterior e, em seguida, a amostragem dos dados pela probabilidade. Em particular, a função de distribuição acima é a preditiva posterior de um dado geométrico com Beta anterior ao parâmetro . p
jII 30/12/2015