Introducción.

En esta publicación se muestra como se implementa de forma computacional el modelo bayesiano Poisson-Gamma y se compara con el resultado analítico.

Modelo Poisson-Gamma.

Se tiene una muestra aleatoria \(X_1\), \(X_2\), \(\ldots\), \(X_n\) con función de distribución de probabilidad Poisson con parámetro \(\lambda\). La función de verosimilitud para este modelo está dada por: \[ \begin{aligned} L\left(\lambda|Datos\right)&=\prod^n_{i=1}\frac{\lambda^{x_i}\exp(-\lambda)}{x_i!}\\ & \propto\lambda^{\sum\limits_{i=1}^nx_i}\exp{(-n\lambda)} \end{aligned} \]

Asumiendo que \(\lambda\) tiene una distribución a priori gamma con parámetros \(\alpha\) y \(\beta\):

\[p\left(\lambda|\alpha,\beta\right)\propto\lambda^{\alpha-1}\exp{\left(-\frac{\lambda}{\beta}\right)} \] Entonces la distribución a posteriori es una gamma con parámetros \(\alpha+\sum\limits^{n}_{i=1}x_i\) y \({\beta}/{(1+\beta n)}\): \[ \begin{aligned} p\left(\lambda|Datos,\alpha,\beta\right)&\propto L\left(\lambda|Datos\right)p\left(\lambda|\alpha,\beta\right)\\ &\propto\lambda^{\alpha+\sum\limits^{n}_{i=1}x_i-1}\exp{\left(-\lambda\frac{(1+\beta n)}{\beta}\right)} \end{aligned} \] El simbolo de proporcionalidad se da porque la expresión presentada en la anterior ecuación no es una función de densidad de probabilidad. Para combertirla en una función de densidad de probabilidad se debe dividir entre \(C\): \[C=\int\limits_{0}^{\infty}\lambda^{\alpha+\sum\limits^{n}_{i=1}x_i-1}\exp{\left(-\lambda\frac{(1+\beta n)}{\beta}\right)} d\lambda \]

Un estimador bayesiano para el parámetro \(\lambda\) está dado por la media de la distribución a posteriori:

\[ \begin{aligned} E\left[\lambda|Datos,\alpha,\beta\right]&=\frac{\int\limits_{0}^{\infty}\lambda^{\alpha+\sum\limits^{n}_{i=1}x_i}\exp{\left(-\lambda\frac{(1+\beta n)}{\beta}\right)} d\lambda }{C}\\ &=\frac{\left(\alpha+\sum\limits^{n}_{i=1}x_i\right)\beta}{(1+\beta n)} \end{aligned} \] A continuación, este estimador se encuentra de manera computacional y se compara con el resultado analítico.

Ilustración

Para ilustrar el modelo bayesiano Poisson-Gamma se asumiran tres escenarios que se distinguen por la cantidad de datos que se tienen:

# Datos Escenario I (n=3)
x <- c(5,4,2)
#Datos Escenario II (n=5)
x1 <- c(5,4,2,1,3)
#Datos Escenario III (n=10)
x2 <- c(5,4,2,1,3,0,2,1,2,2)

En todos estos escenarios se asume \(\alpha=2\) y \(\beta=1\).

Escenario 1

Para obtener el resultado mediante simulación el valor de \(C\) debe ser calculado:

#Numero de simulaciones
m <- 100000
lambdas <- rgamma(n=m, shape=2, scale=1)
aux1 <- function(lambda)prod(dpois(x=x, lambda=lambda))
aux1<-Vectorize(aux1)
C <- mean(aux1(lambdas))

En este caso las curvas de la distribución a priori y a posteriori están dadas por:

posterior <- function(lambda) 
  aux1(lambda)* dgamma(x=lambda, shape=2, scale=1)/C 

prior<-function(lambda){
  dgamma(lambda, shape=2, scale=1)
}
curve(posterior(x), from = 0 , to=9,
      ylab="Density", main="Prior and Posterior Distributions", 
      xlab=expression(lambda), col="blue", las=1,lty=2)
curve(prior(x), from=0, to=9,
      ylab="Density",  
      xlab=expression(lambda), col="tomato", las=1,lty=1,add=TRUE)
points(x=x,y=c(0,0,0),cex=2,col="darkorchid",pch=19)
legend("topright", legend=c("Prior", "Posterior"),
       col=c("red", "blue"), lty=1:2, cex=1.2,border=white)

Los puntos de color morado que aparecen en la figura corresponden a los valores observados en la muestra.

aux2 <- lambdas*aux1(lambdas)
num <- mean(aux2)
num/C

## [1] 3.250967

Este resultado es muy similar al obtenido aplicando la fórmula para la media de la a posteriori que es igual a 3.25.

Escenario 2

En este caso las curvas de la distribución a priori y a posteriori están dadas por:

El estimador bayesiano del parámetro\(\lambda\) es:

aux2 <- lambdas*aux1(lambdas)
num <- mean(aux2)
num/C

## [1] 2.837397

Este resultado es muy similar al obtenido aplicando la fórmula para la media de la a posteriori que es igual a 2.833333.

Escenario 3

Las curvas de la distribución a priori y a posteriori están dadas por:

Comparando esta figura con las anteriores, se observa que a medida que el Tamaño de muestra crece hay una menor variabilidad en la distribución a posteriori. El estimador de \(\lambda\) en este caso es:

aux2 <- lambdas*aux1(lambdas)
num <- mean(aux2)
num/C

## [1] 2.181519

Este resultado es muy similar al obtenido aplicando la fórmula para la media de la a posteriori que es igual a 2.181818.

Modelo Poisson-Gamma