Tengo la siguiente distribución discreta, donde $\alpha,\beta$ son constantes conocidas:
$$ p(x;\alpha,\beta) = \frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)} \;\;\;\;\text{for } x = 0,1,2,\dots $$
¿Cuáles son los métodos más eficaces para tomar muestras de esta distribución?
Dado que $$\frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}=\dfrac{\alpha}{\alpha+\beta+x}\dfrac{\beta+x-1}{\alpha+\beta+x-1}\cdots\dfrac{\beta}{\alpha+\beta}$$ disminuye con $x$ sugiero generar una variante uniforme $u\sim\mathcal{U}(0,1)$ y calculando las sumas acumuladas $$S_k=\sum_{x=0}^k \frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}$$ hasta $$S_k>u$$ La realización es entonces igual a la correspondiente $k$ . Desde $$\eqalign{R_x&=\frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}\\&=\dfrac{\alpha}{\alpha+\beta+x}\dfrac{\beta+x-1}{\alpha+\beta+x-1}\cdots\dfrac{\beta}{\alpha+\beta}\\&=\frac{\alpha+\beta+x-1}{\alpha+\beta+x}\frac{\beta+x-1}{\alpha+\beta+x-1}R_{x-1}\\&=\frac{\beta+x-1}{\alpha+\beta+x}R_{x-1}}$$ y $$S_k=S_k-1+R_k$$ el cálculo puede evitar por completo el uso de funciones Gamma.
(+1) Utilizar $S_k = 1-\frac{\Gamma (a+b) \Gamma (b+k+1)}{\Gamma (b) \Gamma (a+b+k+1)}$ agilizará enormemente el trabajo.
Con respecto a la edición: sospecho que la explotación de las funciones Gamma será útil para resolver el problema de $k$ en términos de $u$ , $\alpha$ y $\beta$ . Por ejemplo, se podría encontrar una aproximación inicial a $u$ utilizando la fórmula de Stirling en la evaluación de $\Gamma(b+k+1)$ y $\Gamma(a+b+k+1)$ y luego pulirlo con algunos pasos de Newton-Raphson. Esos necesitan la evaluación de log Gamma y su derivada. Por supuesto, si $\alpha$ y $\beta$ son ambas integrales entonces la solución es una raíz de un polinomio--pero incluso entonces usar Gamma puede ser el camino a seguir.
¡Gran respuesta! Acepté la respuesta dada por SL porque me llamó la atención un punto clave (que no formaba parte de la pregunta original), que muestrear a partir de la predicción posterior es equivalente a muestrear el parámetro a partir de la posterior, y luego muestrear los datos a partir de la verosimilitud. En particular, la función de distribución anterior es la predictiva posterior de un dato geométrico con una prioridad Beta sobre el parámetro $p$ .
Se trata de un Distribución binomial negativa Beta con el parámetro $r=1$ en tu caso, utilizando la notación de Wikipedia. También nombró Beta-Pascal distribución cuando $r$ es un número entero. Como ha señalado en un comentario, se trata de una distribución predictiva en el modelo binomial negativo bayesiano con una Beta conjugada a priori sobre la probabilidad de éxito.
Por lo tanto, se puede muestrear $\text{Beta}(\alpha,\beta)$ variable $u$ y luego muestrear una variable binomial negativa $\text{NB}(r,u)$ (con $r=1$ en su caso, es decir, una distribución geométrica).
Esta distribución se implementa en el paquete R brr . El muestreador tiene nombre rbeta_nbinom el pmf tiene nombre dbeta_nbinom etc. Las notaciones son $a=r$ , $c=\alpha$ , $d=\beta$ . Compruébalo:
> 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] TRUESi se observa el código, se puede ver que en realidad llama a la función ghyper (hipergeométrica generalizada) de distribuciones de la SuppDists paquete:
brr::rbeta_nbinom
function(n, a, c, d){
rghyper(n, -d, -a, c-1)
}Ineed, la distribución BNB se conoce como una tipo IV distribución hipergeométrica generalizada. Véase la ayuda de ghyper en el SuppDists paquete. Creo que esto también se puede encontrar en el libro de Johnson & al. Distribuciones discretas univariantes .
Esta respuesta está muy bien, pero sería aún mejor si demostraras que la densidad OP publicada es la misma que la densidad binomial negativa.
I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.
