C.2 Métodos de Monte Carlo vía cadenas de Markov

C.2.1 El muestreador de Gibbs

Tal como lo afirma Peña (2002), este procedimiento es apropiado para obtener muestras de una distribución conjunta cuando es fácil muestrear de las distribuciones condicionadas. El algoritmo se implementa asumiendo que \(\boldsymbol \theta_i=(\theta^{(1)}_i, . . . , \theta^{(d)}_i)\) representa a los valores actuales de \(\boldsymbol \theta\). Entonces \(\boldsymbol \theta_{i+1}\) se obtiene así:

  • Generar \(\theta^{(1)}_{i+1}\) de \(p(\theta^{(1)} \mid \theta^{(2)}_i, \ldots,\theta^{(d)}_i,x)\)
  • Generar \(\theta^{(2)}_{i+1}\) de \(p(\theta^{(2)} \mid \theta^{(1)}_{i+1}, \theta^{(3)}_i, \ldots , \theta^{(d)}_i, x)\)
  • \(\ldots\)
  • Generar \(\theta^{(d)}_{i+1}\) de \(p(\theta^{(d)} \mid \theta^{(1)}_{i+1}, \theta^{(2)}_{i+1}, \ldots , \theta^{(d-1)}_{i+1} , x)\)

La idea de este esquema es renovar cada componente por medio de la simulación de la correspondiente distribución condicional. Una vez que la cadena converge, se tiene que los valores de \(\boldsymbol \theta\) corresponden a observaciones de la distribución requerida, \(p(\boldsymbol \theta\mid x)\). Sin embargo, en general, no se garantiza una muestra variables aleatorias totalmente independientes provenientes de la distribución \(p(\theta \mid x)\), dado que el esquema del muestreador de Gibbs usa el valor actual para construir el siguiente valor; por ende, la secuencia de valores que se obtiene estará correlacionada.

Ejemplo C.4 Se puede implementar el muestreador de Gibbs para generar una secuencia de observaciones con densidad conjunta

\[\begin{equation*} (x,y) \sim N_2 \Bigl(0, \begin{pmatrix} \rho & 0 \\ 0 & \rho \end{pmatrix} \Bigl) \end{equation*}\]

Teniendo en cuenta que la media de ambas variables es cero y su varianza uno, entonces la covarianza entre ambas variables será \(\rho\) (Robert y Casella 2009). Por ende, partiendo de valores iniciales \((x_t, y_t)\), el algoritmo se centra en actualizar las distribuciones condicionales según el resultado A.20.

\[\begin{align*} x_{t+1}\mid y_t & \sim N(\rho y_t, 1-\rho^2)\\ y_{t+1}\mid x_{t+1} & \sim N(\rho x_{t+1}, 1-\rho^2) \end{align*}\] 
bivariate.gibbs <- function (n, rho, x, y) {
  mat <- matrix(ncol = 2, nrow = n)
  mat[1, ] <- c(x, y)
  for (i in 2:n){
    x <- rnorm(1, rho * y, sqrt(1 - rho^2))
    y <- rnorm(1, rho * x, sqrt(1 - rho^2))
    mat[i, ] <- c(x, y)
  }
  mat<-as.data.frame(mat)
  return(mat)
}

biv <- bivariate.gibbs(n=2000, rho=0.5, x= 0, y = 0)
colMeans(biv)
##          V1          V2 
## -0.02273293 -0.01762626
var(biv)
V1 V2
V1 0.9702333 0.4709515
V2 0.4709515 0.9728540 
cor(biv)
V1 V2
V1 1.000000 0.484746
V2 0.484746 1.000000 
plot(biv)
Generación de valores para una distribución normal bivariada.

Figura C.7: Generación de valores para una distribución normal bivariada.

Ejemplo C.5 Un problema común es el de descartar los primeros valores, puesto que el algoritmo puede demorar en obtener convergencia; esto se puede resolver en forma empírica utilizando las medias y varianzas acumuladas y graficándolas se puede tomar una decisión acerca del valor óptimo en el que la cadena converge.

Con el siguiente código computacional, es posible corroborar que un punto de corte óptimo desde el cual se consideraría que las cadenas simuladas anteriormente es a partir de la iteración 600. 
g.diag <- function(sample){
  n <- length(sample) 
  res <- matrix(nrow=2, ncol=n)
  for(i in 1:n){
    res[1, i] <- mean(sample[1 : i])
    res[2, i] <- var(sample[1 : i])
  }
  return(res)
}

m1 <- g.diag(biv[, 1])
m2 <- g.diag(biv[, 2])

par(mfcol = c(1, 2))
plot(m1[1, ], type = 'l', ylim=c(-0.6, 0.6), col=4)
lines(m2[1, ], lty = 2, col = 2)
title("Diagnóstico - Media acumulada")

plot(m1[2, ], type = 'l', ylim = c(0.5, 1.5), col=4)
lines(m2[2, ], lty = 2, col = 2)
title("Diagnóstico - Varianza acumulada")
Convergencia de la media y varianza usando el muestreador de Gibbs.

Figura C.8: Convergencia de la media y varianza usando el muestreador de Gibbs.

El muestrador de Gibbs también funciona en una “segunda fase”, cuando queremos seleccionar una muestra de \(f(\theta\mid x)\), es decir, la distribución de los parámetros dada la información observada \(x\).

Ejemplo C.6 Suponga que \(y\) tiene distribución \(N(\mu,\sigma^2=1/\phi)\) y queremos obtener una muestra de la distribución posterior del vector aleatorio \(\boldsymbol \theta=(\mu,1/\phi)\). Para este caso supongamos que conocemos las distribuciones previas; para la media \(\mu\) se asume una distribución uniforme y para la varianza \(\phi\) una distribución Gamma con parámetros \(a\) y \(b\). La distribución posterior de \((\mu, \phi)\) satisface:

\[\begin{equation} p(\mu, \phi \mid y) \propto (\phi)^{n/2} \exp\left\{-\phi \frac{\sum_{j=1}^n(y_j-\mu)^2}{2}\right\}(\phi)^{a-1}exp(-b/\phi) \end{equation}\]

En donde la primera parte después del signo de proporcionalidad, corresponde a la verosimilitud de la información observada y la segunda parte corresponde a la distribución posterior de \(\phi\); la distribución posterior de \(\mu\) no aparece pues es una constante. Por tanto, ésta se puede escribir como:

\[\begin{equation*} p(\mu, \phi \mid y)\propto(\phi)^{n/2+a-1}exp\left\{-\phi\Bigl(\frac{\sum_{j=1}^n(y_j-\mu)^2}{2}+b\Bigl)\right\} \end{equation*}\]

Acudiendo al resultado A.10, la distribución condicional de la varianza \(\sigma^2\) dado \((\mu, y)\) es Gamma-inversa con parámetros \(a+n/2\) y \(\sum_{j=1}^n(y_j-\mu)^2/2+b\). Por tanto,

\[\begin{equation} \tag{C.1} \sigma^2\mid\mu,x\sim Gamma-inversa\biggl(\theta+n/2,\sum_{j=1}^n(y_j-\mu)^2/2+b\biggl) \end{equation}\]

Análogamente, la distribución de \(\mu\) dado \((\sigma^2, y)\) es normal con media \(\bar{y}\) y varianza \(\sigma^2/n\), es decir,

\[\begin{equation} \tag{C.2} \mu\mid\sigma^2,y\sim N(\bar{y},\sigma^2/n) \end{equation}\]

Para implementar el muestreador de Gibbs con estas distribuciones, primero se deben escoger valores apropiados para \(a\) y \(b\), con el propósito de representar correctamente la distribución previa, y luego

  • Defninir un valor inicial para la media y la varianza, \((\mu_0, \sigma^2_0)\).
  • Generar \((\mu_{i+1}, \sigma_{i+1}^2)\) simulando \(\mu_{i+1}\) de (C.1) y luego \(\sigma^2_{i+1}\) de (C.2).
  • Iterar para obtener \((\mu_0, \sigma^2_0), (\mu_1, \sigma^2_1), (\mu_2, \sigma^2_2),\cdots,\).
  • Suponiendo que el algoritmo converge después de \(m\) iteraciones, descartar los \(m\) primeros valores.
Entonces \((\mu_{m+1}, \sigma^2_{m+1}), (\mu_{m+2}, \sigma^2_{m+2}),\cdots,\) es una muestra (correlacionada) de \(p(\mu, \sigma^2\mid x)\).

La siguiente función en R implementa el muestreador de Gibbs para el anterior ejemplo.

library(invgamma)

normal2 <- function(datos, a, b, nsim, inicial){
  n <- length(datos)
  xbar <- mean(datos)
  mu.now <- inicial[1]
  var.now <- inicial[2]
  dummy <- matrix(ncol = 2, nrow = nsim)
  dummy[1, 1] <- mu.now
  dummy[1, 2] <- var.now
  
  for (i in 2 : nsim){
    alp <- a + (n/2)
    bet <- b + (sum((datos - mu.now)^2)/2)
    var.next <- rinvgamma(1, shape = alp, rate = bet)
    mu.next <- rnorm(1, xbar, sqrt(var.now/n))
    dummy[i, 1] <- mu.next
    dummy[i, 2] <- var.next
    mu.now <- mu.next
    var.now <- var.next
  }
  return(dummy)
}

datos <- rnorm(100, 5, 2)
mc1.vals <- normal2(datos, a = 2, b = 5, 
                    nsim = 1000, inicial = c(2, 2))
mc1.vals <- mc1.vals[101: 1000, ]
colMeans(mc1.vals)
## [1] 5.257333 4.107664
par(mfcol = c(1, 2))
plot(mc1.vals[, 1], type = 'l', ylab = 'mu')
plot(mc1.vals[, 2], type = 'l', ylab = 'sigma^2')
Cadenas generadas desde el muestreador de Gibbs.

Figura C.9: Cadenas generadas desde el muestreador de Gibbs. 

par(mfcol = c(1, 2))
hist(mc1.vals[, 1], prob = T, xlab='mu', main = "")
lines(density(mc1.vals[, 1], kernel='gaussian'))
hist(mc1.vals[, 2], prob = T, xlab='sigma^2', main = "")
lines(density(mc1.vals[, 2], kernel='gaussian'))
Densidades posteriores generadas con el muestreador de Gibbs.

Figura C.10: Densidades posteriores generadas con el muestreador de Gibbs.

C.2.2 El algoritmo de Metrópolis-Hastings

Este algoritmo se basa en proponer un nuevo punto de acuerdo a una función de densidad adecuada y aceptar este nuevo valor propuesto con una probabilidad que depende del punto actual, del nuevo punto y de la densidad de la cual fue propuesto el nuevo punto.

Suponga que deseamos simular valores de una distribución multivariada \(p(\theta \mid y)\). Sea la función de densidad propuesta \(q(\theta, \theta')\), una función de densidad de probabilidad arbitraria que describe la probabilidad de aceptación de \(\theta'\) a partir de la posición actual de \(\theta\). El algoritmo de Metropolis-Hastings está dado por los siguientes pasos:

  • Siendo el valor actual \(\theta_i\), genere un valor candidato \(\theta'\) obtenido como una observación de la densidad \(q(\theta_i, \theta')\).
  • Calcule \[\begin{equation*} T(\theta_i, \theta') = \begin{cases} \min \left(1, \frac{p(\theta' \mid y)q(\theta', \theta_i)}{p(\theta_i \mid y)q(\theta_i, \theta')} \right), & \text{ si } p(\theta_i \mid y)q(\theta_i, \theta') > 0,\\ 1, & \text{ si } p(\theta_i \mid y)q(\theta_i, \theta') = 0 \end{cases} \end{equation*}\]
  • Acepte el nuevo valor y actualícelo a \(\theta_{i+1}=\theta'\) con probabilidad \(T(\theta_i, \theta')\). De otra forma, rechazar el valor candidato y defina \(\theta_{i+1}=\theta_i\). Repita el paso anterior para obtener la secuencia \(\theta_0,\theta_1,...,\) donde \(\theta_0\) denota un valor arbitrario de arranque. Descarte los primeros \(m\) valores obtenidos.

Siguiendo el anterior algoritmo, entonces se tiene que \(\theta_{m+1}, \theta_{m+2}, \ldots\) es una secuencia (correlacionada) de la distribución requerida. En principio, puede ser usada cualquier densidad \(q\), pero si ésta es escogida ingenuamente, la eficiencia de la cadena puede ser muy pobre. La relación más importante entre el muestreados de Gibbs y el algorítmo de Metropolis-Hastings, está dada como un teorema en el libro de Robert y Casella (2009pág. 296).

Resultado A.7 El muestreador de Gibbs es equivalente al algoritmo de Metropolis-Hastings, con la probabilidad de aceptación igual a uno para todos los puntos propuestos.


Lo anterior implica que la convergencia para ambos métodos no es la misma. Para cerrar la sección de cadenas de Markov vía Monte Carlo, es importante hacernos la siguiente pregunta: ¿Son independientes las muestras simuladas? En principio no se puede hablar de independencia, pues es claro que la observación \(\{i+1\}\) depende de la observación \(\{i\}\). Dado que las observaciones resultantes se encuentran en estricto orden de medición, podríamos utilizar algunos criterios como la función de auto-correlación (ACF) y la función de auto-correlación parcial (PACF), para conocer sobre la correlación entre observaciones.

Siguiendo con el ejemplo C.6 del apartado de Gibbs, se ha escogido usar como como distribuciones propuestas \(q\) para la media y para la varianza, densidades normales centradas en el actual parámetro, ambas con varianza igual a uno. Dadas las distribuciones propuestas, algunos valores de la varianza pueden ser negativos; aunque este no es un problema porque la distribución posterior le asignará el valor cero, por tanto este valor será rechazado con un probabilidad de uno.

library(invgamma)

met.hast <- function(datos, a, b, iter, ini){
  mu0 <- ini[1] 
  var0 <- ini[2]
  resul <- matrix(ncol = 2, nrow = iter)
  resul[1, 1] <- mu0
  resul[1, 2] <- var0
  for (i in 2 : iter){
    mu.prop <- rnorm(1, mu0, 1)
    var.prop <- rnorm(1, var0, 1)
    if (var.prop <= 0){ T.val <- 0 }
    else{
      p1 <- prod(dnorm(datos, mu.prop, sqrt(var.prop))) *
        dinvgamma(var.prop, shape = a, rate = b)
      q1 <- dnorm(mu0, mu.prop, 1) *
        dnorm(var0, var.prop, 1)
      p2 <- prod(dnorm(datos, mu0, sqrt(var0))) * 
        dinvgamma(var0, shape = a, rate = b)
      q2 <- dnorm(mu.prop, mu0, 1) * 
        dnorm(var.prop, var0, 1)
      T.val <- min(1, (p1 * q1)/(p2 * q2))
    }
    u <- runif(1) 
    if (u <= T.val){
      resul[i, 1] <- mu.prop
      resul[i, 2] <- var.prop
      } 
    else{
      resul[i, 1] <- mu0
      resul[i, 2] <- var0
      }
    mu0 <- resul[i, 1]
    var0 <- resul[i, 2]
  }
  return(resul)
}

datos <- rnorm(100, 5, 2)
mc2 <- met.hast(datos, a = 2, b = 5, 
                iter = 1000, ini = c(2, 2))
colMeans(mc2)
## [1] 5.061780 4.744852
par(mfrow=c(2,2))
pacf(mc2[, 1], 100)
pacf(mc2[, 2], 100)
acf(mc2[, 1], 100)
acf(mc2[, 2], 100)
Autocorrelación y autocorrelación parcial para las cadenas simuladas del algoritmo MH.

Figura C.11: Autocorrelación y autocorrelación parcial para las cadenas simuladas del algoritmo MH.

C.2.3 Buenas prácticas en la aplicación de métodos MCMC

Dado que una gran parte de la inferencia bayesiana está ligada a la programación e implementación de los métodos MCMC para realizar inferencias posteriores de los parámetros de interés, se sugiere seguir el razonamiento y recomendaciones de A. Gelman y Shirley (2010), que puede ser resumido en los siguientes ítemes para cada parámetro de interés:

  1. Simulación de tres o más cadenas de forma paralela. Los valores iniciales de cada cadena deben estar dispersos entre sí.
  2. Comprobación de la convergencia de las cadenas mediante el descarte de la primera mitad de los valores generados en las cadenas. Esta etapa se conoce como burning stage.
  3. Una vez que las cadenas converjan, mezclar los tres conjuntos de valores generados por las cadenas. Esto garantiza, en primera instancia, que las cadenas no estén auto-correlacionadas.
  4. Además de realizar esta mezcla, descartar valores intermedios mediante un muestreo sistemático. Esta etapa se conoce como thining stage. Al final se recomienda almacenar una cantidad elevada de valores simulados.
  5. Calibrar el algoritmo si la convergencia de las cadenas no se presenta rápidamente.
    • Para los algoritmos de Metropolis-Hastings, escoger una distribución de salto acorde con la distribución de la cual se desea simular. Por ejemplo, Cepeda y Gamerman (2001) presentan dos distribuciones de salto para el problema de la modelación de la varianza (cada una de las propuestas presenta tasas de aceptación diferentes).
  6. Comparación y contraste de los resultados con modelos simples que permitan examinar posibles discrepancias y corregir errores de programación.

En términos de inferencia bayesiana, se tienen dos tipos de procesos: el primero y más común, que trata de realizar inferencias acerca de un vector de parámetros de interés \(\boldsymbol \theta\); el segundo trata con los momentos del parámetro, por ejemplo su esperanza. Nótese que el primer proceso se presenta con seguridad en ejercicios empíricos simulados; sin embargo, el segundo se presenta en los ejercicios prácticos con datos reales, en donde se quiere contrastar alguna hipótesis.

Las anteriores dos opciones tienen tratamientos muy diferentes en términos de la cantidad de simulaciones requeridas. Por ejemplo, si el objetivo es inferir acerca de \(\boldsymbol \theta\), para conocer su comportamiento estructural, basta con realizar una simulación que genere una cantidad mediana de valores y que se resumen en un promedio y una desviación estándar. Por otro lado, si el objetivo es inferir acerca de \(E(\boldsymbol \theta)\), se requieren muchas más simulaciones para obtener una buena precisión. Siguiendo a A. Gelman y Shirley (2010), una vez terminado el proceso de burning y thining, se sugiere que se dividan los valores simulados en las cadenas paralelas y se formen \(k\) grupos; de esta forma, una estimación de \(E(\boldsymbol \theta)\) será la gran media de las medias muestrales de cada grupo y el error estándar será su desviación estándar dividida por \(\sqrt{k}\).

Referencias

Cepeda, E., y D. Gamerman. 2001. «Bayesian Modeling of Variance Heterogeneity in Normal Regression Models». Brazilian Journal of Probability and Statistics 14: 207-21.
Gelman, A., y K Shirley. 2010. «Handbook of Markov Chain Monte Carlo». En. CRC.
Peña, D. 2002. Análisis de datos multivariantes. McGraw-Hill.
———. 2009. Introducing Monte Carlo Methods with R. Springer.