Muestreador de Gibbs

1 Introducción

En la mayoría de los casos no es posible obtener directamente muestras la distribución posterior.

2 Motivación: modelo Normal con una distribución previa semi-conjugada

Considere especificar Su estado de información previo acerca de \(\theta\) de manera independiente de \(\sigma^2\) de forma que: \[ p(\theta,\sigma^2) = p(\theta)\,p(\sigma^2)\,. \] Esta formulación es mas flexible porque no hay una restricción de dependencia a priori entre \(\theta\) y \(\sigma^2\).

Asumiendo que las observaciones dadas en \(\boldsymbol{y} = (y_1,\ldots,y_n)\) son intercambiables, bajo el modelo Normal con una distribución previa semi-conjugada se tiene que:

• Verosimilitud (distribución muestral conjunta): \[ y_i\mid\theta,\sigma^2 \stackrel{\text{iid}}{\sim} \textsf{N}(\theta,\sigma^2)\,,\qquad i=1,\ldots,n. \]

• Previa: \[\begin{align*} \theta &\sim \textsf{N}(\mu_0, \tau^2_0) \\ \sigma^2 &\sim \textsf{GI}\left(\tfrac{\nu_0}{2},\tfrac{\nu_0\,\sigma^2_0}{2}\right) \end{align*}\]

• Hiper-parámetros: \(\mu_0\), \(\tau^2_0\), \(\nu_0\), y \(\sigma^2_0\).

En el caso de la previa conjugada (donde \(\tau_0^2\) es proporcional a \(\sigma^2\)) se tiene que la densidad marginal posterior de \(\sigma^2\) es Gamma Inversa. En este caso, \(\sigma^2\mid \boldsymbol{y}\) no sigue una distribución estándar conocida de la cual se pueda obtener muestras fácilmente.

2.1 Distribuciones condicionales completas

El muestreador de Gibbs (Gibs sampler) es un algoritmo iterativo que permite obtener muestras (correlacionadas) de la distribución posterior por medio de las distribuciones condicionales completas.

Bajo esta especificación del modelo Normal, se demuestra que:

• La distribución condicional completa de \(\theta\) es \(\theta\mid\sigma^2,\boldsymbol{y}\sim\textsf{N}(\mu_n,\tau^2_n)\), donde \[ \mu_n = \frac{\frac{1}{\tau^2_0}\mu_0 + \frac{n}{\sigma^2}\bar{y}}{\frac{1}{\tau^2_0} + \frac{n}{\sigma^2}} \qquad\text{y}\qquad\tau^2_n=\frac{1}{\frac{1}{\tau^2_0} + \frac{n}{\sigma^2}}\,. \]
• La distribución condicional completa de \(\sigma^2\) es \(\sigma^2\mid\theta,\boldsymbol{y}\sim\textsf{GI}\left(\tfrac{\nu_n}{2},\tfrac{\nu_n\,\sigma^2_n}{2}\right)\), donde \[ \nu_n = \nu_0+n\qquad\text{y}\qquad\sigma^2_n = \frac{1}{\nu_n}\left( \nu_0\sigma^2_0 + ns^2(\theta) \right)\,. \] con \(s^2(\theta) = \tfrac{1}{n}\sum_{i=1}^n (y_i-\theta)^2 = (n-1)s^2 + n(\bar{y}-\theta)^2\), el cual corresponde al estimardor insesgado de \(\sigma^2\) si \(\theta\) fuera conocido. Es más eficiente calcular la segunda expresión.

3 Muetreador de Gibbs

Dado un estado actual de los parámetros del modelo \(\boldsymbol{\theta}^{(b)} = \left(\theta^{(b)}, (\sigma^2)^{(b)}\right)\), se genera un nuevo estado \(\boldsymbol{\theta}^{(b+1)}\) como sigue:

1. Muestrar \(\theta^{(b+1)}\sim p(\theta\mid(\sigma^2)^{(b)}, \boldsymbol{y})\).
2. Muestrar \((\sigma^2)^{(b+1)}\sim p(\sigma^2\mid\theta^{(b+1)}, \boldsymbol{y})\).
3. Establecer \(\boldsymbol{\theta}^{(b+1)} = \left(\theta^{(b+1)}, (\sigma^2)^{(b+1)}\right)\).
4. Repetir los pasos 1. a 3. hasta convergencia.

Este algoritmo se denomina muestreador de Gibbs y genera una secuencia dependiente de parámetros \(\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(B)}\) de la distribución posterior conjunta \(p(\theta,\sigma^2\mid \boldsymbol{y})\).

Como punto de partida, solo es necesario proporcionar un valor inicial para \(\sigma^2\). Usualmente este valor se muestra de la distribución previa correspondiente, esto es, \((\sigma^2)^{(0)}\sim\textsf{IG}\left(\tfrac{\nu_0}{2},\tfrac{\nu_0\,\sigma^2_0}{2}\right)\).

3.1 Algoritmo general

Dado un estado actual de los parámetros del modelo \(\boldsymbol{\theta}^{(b)} = \left(\theta_1^{(b)},\ldots,\theta_k^{(b)}\right)\), se genera un nuevo estado \(\boldsymbol{\theta}^{(b+1)}\) como sigue:

1. Muestrear \(\theta_1^{(b+1)}\sim p\left(\theta_1\mid\theta_2^{(b)},\theta_3^{(b)},\ldots,\theta_k^{(b)}\right)\).
2. Muestrear \(\theta_2^{(b+1)}\sim p\left(\theta_2\mid\theta_1^{(b+1)},\theta_3^{(b)},\ldots,\theta_k^{(b)}\right)\).
3. Muestrear \(\theta_3^{(b+1)}\sim p\left(\theta_3\mid\theta_1^{(b+1)},\theta_2^{(b+1)},\ldots,\theta_k^{(b)}\right)\).
4. …
5. Muestrear \(\theta_k^{(b+1)}\sim p\left(\theta_3\mid\theta_1^{(b+1)},\theta_2^{(b+1)},\ldots,\theta_{k-1}^{(b+1)}\right)\).
6. Establecer \(\boldsymbol{\theta}^{(b+1)} = \left(\theta_1^{(b+1)},\ldots,\theta_k^{(b+1)}\right)\).
7. Repetir los pasos 1. a 6. hasta convergencia.

Este algoritmo se denomina muestreador de Gibbs y genera una secuencia dependiente de parámetros \(\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(B)}\) de la distribución posterior conjunta \(p(\theta_1,\ldots,\theta_k\mid \boldsymbol{y})\).

Observe que \(\boldsymbol{\theta}^{(b+1)}\) depende únicamente de \(\boldsymbol{\theta}^{(b)}\) lo cual sugiere que \(\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(B)}\) define una cadena de Markov.

Se puede demostrar que la cadena converge a la distribución objetivo (distribución posterior) a medida que \(b\rightarrow B\) sin importar el punto de partida \(\boldsymbol{\theta}^{(0)}\) (algunos puntos de partida pueden ser mejores que otros). Por lo tanto es posible usar la secuencia \(\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(B)}\) tal y como se tratara de una muestra aleatoria de valores de \(\boldsymbol{\theta}=(\theta_1,\ldots,\theta_k)\) provenientes de la distribución posterior \(p(\theta_1,\ldots,\theta_k\mid y)\). Estos es, \[ \frac{1}{B} \sum_{b=1}^{B} g(\boldsymbol{\theta}) \longrightarrow \textsf{E}[g(\boldsymbol{\theta}) \mid \boldsymbol{y}]=\int_{\Theta} g(\boldsymbol{\theta}) p(\boldsymbol{\theta} \mid \boldsymbol{y})\qquad\text{cuando}\qquad b\rightarrow \infty \]

El punto de partida \(\boldsymbol{\theta}^{(0)}\) es arbitrario y usualmente se muestrea de la distribución previa.

El muestreador de Gibbs hace parte de un conjunto de técnicas de aproximación denominadas Markov Chain Monte Carlo (MCMC).

Los métodos de MCMC son técnicas de aproximación numérica, no son modelos, no generan más información además de la contenida en \(\boldsymbol{y}\) y \(p(\boldsymbol{\theta})\), simplemente son métodos numéricos.

3.2 Diagnosticos de convergencia (estacionaridad)

Para que las aproximaciones a las cantidad posteriores sean precisas, necesitamos que la distribución empírica de \(\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(B)}\) esté lo suficientemente cerca de \(p(\boldsymbol{\theta}\mid y)\).

El muestreador de Gibbs siempre produce muestras que eventualmente van a converger a la distribución objetivo, pero en algunos casos la convergencia puede ser lenta debido a la autocorrelación de los parámetros.

Es usual pensar en la secuencia \(\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(B)}\) como la trayectoria de una partícula \(\boldsymbol{\theta}\) moviéndose a lo largo del espacio de parámetros \(\Theta\).

En general, nunca es posible estar absolutamente seguro que la cadena a convergido. Pero por lo menos podemos saber si no lo ha hecho.

3.2.1 Traceplots

Permiten chequear estacionaridad (muestras de una parte de la cadena tienen una distribución similar a muestras de otra parte de la cadena).

También permiten ver si la cadena está mezclando bien o no (en simulación de Monte Carlo la partícula se mueve libremente a cualquier región del espacio de parámetros, lo cual significa cero autocorrelación).

Pare resolver problemas de estacionaridad se recomienda correr la cadena con más iteraciones.

3.2.2 Autocorrelación

La función de autocorrelación está dada por \[ \operatorname{acf}_{t}(\theta)=\frac{\frac{1}{B-t} \sum_{b=1}^{B-t}\left(\theta^{(b)}-\bar{\theta}\right)\left(\theta^{(b+t)}-\bar{\theta}\right)}{\frac{1}{B-1} \sum_{b=1}^{B}\left(\theta^{(b)}-\bar{\theta}\right)^{2}}\qquad\text{donde}\qquad\bar{\theta} = \frac{1}{B}\sum_{b=1}^B \theta^{(b)}\,. \] Entre mayor sea la autocorrelación, se necesitan más muestras para obtener la precisión deseada.

El tamaño efectivo de la muestra para una secuencia de muestras de \(\boldsymbol{\theta}\) obtenidos usando MCMC.

4 Ejemplo

En 1981, los biólogos W. L. Grogan y W. W. Wirth descubrieron en las selvas de Brasil dos nuevas variedades de un diminuto insecto picador llamado mosquito (midge). Llamaron a un tipo de mosquito Apf y al otro mosquito Af. Los biólogos descubrieron que el mosquito Apf es portador de una enfermedad debilitante que causa inflamación del cerebro cuando un humano está mordido por un mosquito infectado. Aunque la enfermedad rara vez es fatal, la discapacidad causada por la hinchazón puede ser permanente. La otra forma de mosquito, el Af, es bastante inofensiva y un valioso polinizador. En un esfuerzo por distinguir las dos variedades, los biólogos tomaron medidas en los mosquitos que capturaron. Este es un conjunto de datos valioso para probar métodos de clasificación.

Considere los datos de la longitud del ala en milímetros (\(y\)) de \(n=9\) miembros de la especie Af de mosquitos. A partir de estas nueve mediciones, se quiere hacer inferencia sobre la media poblacional \(\theta\). Otros estudios sugieren que la longitud de las alas suele ser de alrededor de 1.9 mm. Claramente, se tiene que las longitudes deben ser positivas, lo que implica que \(\theta > 0\).

Los datos de Grogan y Wirth se encuentran disponibles en la libreria Flury de R, pero esta libraria no se encuentra disonible para versiones resientes de R.

# datos 
y <- c(1.64,1.70,1.72,1.74,1.82,1.82,1.82,1.90,2.08)

# tamaño de la muestra
n <- length(y)

# estadisticos
mean.y <- mean(y)
var.y  <- var(y)
sum.y  <- sum(y)

4.1 Muestreador de Gibbs para el modelo Normal

# hiperparametros
mu0 <- 1.9 
t20 <- 0.5^2
s20 <- 0.01 
nu0 <- 1

# numero de muestras 
S <- 100000
# matriz para almacenar las muestras
PHI <- matrix(data = NA, nrow = S, ncol = 2)
# mostrar anuncios cada 10% de las iteraciones
ncat <- floor(S/10) 

# ALGORITMO (muestreador de Gibbs)
# 1. inicializar la cadena
#    valor inicial: simular de la previa
#    solo es necesario alguno de los valores
set.seed(1)
phi <- c( rnorm(1, mu0, sqrt(t20)), rgamma(1, nu0/2, nu0*s20/2) )
PHI[1,] <- phi
# 2. simular iterativamente de las distribuciones condicionales completas
set.seed(2)
for(s in 2:S) {
        # 2.1 actualizar el valor de theta
        t2n    <- 1/( 1/t20 + n*phi[2] )      
        mun    <- (mu0/t20 + sum.y*phi[2]) * t2n
        phi[1] <- rnorm( 1, mun, sqrt(t2n) )
        # 2.2 actualizar el valor de sigma^2
        nun    <- nu0+n
        s2n    <- (nu0*s20 + (n-1)*var.y + n*(mean.y-phi[1])^2)/nun
        phi[2] <- rgamma(1, nun/2, nun*s2n/2)
        # 2.3 almacenar
        PHI[s,] <- phi
        # 2.4 progreso
        if (s%%ncat == 0) cat(100*round(s/S, 1), "% completado ... \n", sep = "")
}

## 10% completado ... 
## 20% completado ... 
## 30% completado ... 
## 40% completado ... 
## 50% completado ... 
## 60% completado ... 
## 70% completado ... 
## 80% completado ... 
## 90% completado ... 
## 100% completado ...

# grafico del algoritmo
# este grafico no se acostumbra a hacer en la practica

par(mfrow=c(1,3),mar=c(2.75,3,.5,.5),mgp=c(1.70,.70,0))
m1<-5
plot(PHI[1:m1,],type="l",xlim=range(PHI[1:100,1]), ylim=range(PHI[1:100,2]),
      lty=1,col="gray",xlab=expression(theta),ylab=expression(tilde(sigma)^2))
text(PHI[1:m1,1], PHI[1:m1,2], c(1:m1) )

m1<-15
plot(PHI[1:m1,],type="l",xlim=range(PHI[1:100,1]), ylim=range(PHI[1:100,2]),
      lty=1,col="gray",xlab=expression(theta),ylab=expression(tilde(sigma)^2))
text(PHI[1:m1,1], PHI[1:m1,2], c(1:m1) )

m1<-100
plot(PHI[1:m1,],type="l",xlim=range(PHI[1:100,1]), ylim=range(PHI[1:100,2]),
      lty=1,col="gray",xlab=expression(theta),ylab=expression(tilde(sigma)^2))
text(PHI[1:m1,1], PHI[1:m1,2], c(1:m1) )

4.2 Diagnosticos de convergencia

# log-verosimilitud
LL <- NULL
for (i in 1:S)
  LL[i] <- sum(dnorm(x = y, mean = PHI[i,1], sd = sqrt(1/PHI[i,2]), log = T))

# graficos
par(mfrow=c(3,1), mar=c(3,3,1,1), mgp=c(1.75,.75,0))
# traceplots
plot(PHI[,1], type = "l", xlab="iteration", ylab=expression(theta))
abline(h = mean(PHI[,1]), col = 4, lty = 2, lwd = 2)
plot(1/PHI[,2], type = "l", xlab="iteration", ylab=expression(sigma^2))
abline(h = mean(1/PHI[,2]), col = 4, lty = 2, lwd = 2)
plot(LL, type = "l", xlab="iteration", ylab="lLog-verosimilitud")
abline(h = mean(LL), col = 4, lty = 2, lwd = 2)

# autocorrelacion
par(mfrow = c(1,2), mar=c(3,3,3,1), mgp=c(1.75,.75,0))
acf(PHI[,1], main = expression(theta))
acf(1/PHI[,2], main = expression(sigma^2))
# tamaño efectivo de la muestra
coda::effectiveSize(PHI)

##     var1     var2 
## 100000.0  84767.8

4.3 Distribuciones posterior y marginales

par(mfrow=c(2,2),mar=c(2.75,3,.5,.5),mgp=c(1.70,.70,0))
sseq <- 1:10000
# distribucion conjunta
plot(PHI[sseq,1], PHI[sseq,2], pch=".", xlim=range(PHI[,1]),ylim=c(0,225), xlab=expression(theta), ylab=expression(tilde(sigma)^2))
# distribucion conjunta
plot(PHI[sseq,1], 1/PHI[sseq,2], pch=".", xlim=range(PHI[,1]),ylim=c(0,0.15), xlab=expression(theta), ylab=expression(sigma^2))
# theta
plot(density(PHI[,1], adj=2), xlab=expression(theta), main="", xlim=c(1.55,2.05), ylab=expression(paste(italic("p("), theta,"|",italic(y[1]),"...",italic(y[n]),")",sep="")))
# precision
plot(density(PHI[,2],adj=2), xlab=expression(tilde(sigma)^2),main="", ylab=expression(paste(italic("p("),tilde(sigma)^2,"|",italic(y[1]),"...",italic(y[n]),")",sep=""))) 

4.4 Inferencia

# intervalos de credibilidad
# media
round(quantile(PHI[,1], c(.025, 0.5, 0.975)), 3)

##  2.5%   50% 97.5% 
## 1.710 1.805 1.901

# precision 
round(quantile(PHI[,2], c(.025, 0.5, 0.975)), 3)

##    2.5%     50%   97.5% 
##  18.857  57.758 131.302

# desviancion estandar
round(quantile(1/sqrt(PHI[,2]), c(.025, 0.5, 0.975)), 3)

##  2.5%   50% 97.5% 
## 0.087 0.132 0.230

# coeficiente de variacion
round(quantile((1/sqrt(PHI[,2]))/PHI[,1], c(0.025, 0.5, 0.975)), 3)

##  2.5%   50% 97.5% 
## 0.048 0.073 0.128

# probabilidad posterior
round(mean(PHI[,1] > 1.8), 3)

## [1] 0.547