8.2 Cadenas de Markov Monte Carlo

Markov Chain Monte Carlo (MCMC) es una familia de algoritmos que generan muestras de una distribución objetivo construyendo una cadena de Markov cuya distribución estacionaria es la distribución a posteriori. Las muestras están correlacionadas pero, bajo condiciones de regularidad, convergen a la distribución a posteriori.

8.2.1 ¿Por qué MCMC?

Hasta ahora, usamos la aproximación cuadrática (quap) para estimar la distribución a posteriori. Esto asume que la distribución a posteriori es aproximadamente normal multivariada. Cuando este supuesto falla — distribuciones a posteriori no normales, parámetros discretos, dependencias complejas — necesitamos MCMC.

MCMC funciona para cualquier distribución a posteriori, sin importar su forma, porque no aproxima la distribución a posteriori — la muestrea.

$p(\theta \mid D) = \frac{p(D \mid \theta) \cdot p(\theta)}{p(D)}$

MCMC evita calcular $p(D)$ (la constante de normalización) usando el cociente de probabilidades no normalizadas — igual que el rey Markov.

# Compare quap (normal approximation) with MCMC samples
# for a posterior that is NOT normal

set.seed(42)

# Simulate from a skewed posterior (exponential-like)

# True posterior: theta ~ Gamma(2, 1)

n_samples <- 10000

# MCMC samples (here we just draw from the true distribution)

mcmc_samples <- rgamma(n_samples, shape = 2, rate = 1)

# Quadratic approximation: mean and SD of the samples

mean_q <- mean(mcmc_samples)
sd_q <- sd(mcmc_samples)
quap_samples <- rnorm(n_samples, mean = mean_q, sd = sd_q)

cat("Comparison: MCMC vs Quadratic approximation\n")
cat(" MCMC mean:", round(mean(mcmc_samples), 3), " SD:", round(sd(mcmc_samples), 3), "\n")
cat(" quap mean:", round(mean(quap_samples), 3), " SD:", round(sd(quap_samples), 3), "\n")
cat(" MCMC 95% CI: [", round(quantile(mcmc_samples, 0.025), 3), ",",
round(quantile(mcmc_samples, 0.975), 3), "]\n")
cat(" quap 95% CI: [", round(quantile(quap_samples, 0.025), 3), ",",
round(quantile(quap_samples, 0.975), 3), "]\n")

# Plot

par(mfrow = c(1, 2))
hist(mcmc_samples, breaks = 50, col = "steelblue", border = "white",
main = "MCMC (true posterior)", xlab = "theta", freq = FALSE)
hist(quap_samples, breaks = 50, col = "firebrick", border = "white",
main = "Quadratic approximation", xlab = "theta", freq = FALSE)
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

8.2.2 El algoritmo de Metropolis-Hastings

El algoritmo de Metropolis-Hastings generaliza el algoritmo de Metropolis para permitir distribuciones de propuesta arbitrarias:

Proponer: Dibujar $\theta' \sim q(\theta' \mid \theta)$ de una distribución de propuesta.
Cociente de aceptación: $R = \frac{p(\theta' \mid D) \cdot q(\theta \mid \theta')}{p(\theta \mid D) \cdot q(\theta' \mid \theta)}$
Aceptar/Rechazar: Aceptar con probabilidad $\min(1, R)$ .

# Metropolis-Hastings for a simple Bayesian model
# Model: y ~ Normal(mu, 1), prior: mu ~ Normal(0, 10)
# Posterior: mu ~ Normal(ybar, 1/n) (conjugate)

set.seed(42)
n <- 20
y <- rnorm(n, mean = 3, sd = 1) # true mu = 3
ybar <- mean(y)

# Metropolis-Hastings

n_iter <- 5000
mu <- 0 # starting value
samples <- numeric(n_iter)

for (t in 1:n_iter) {

# Propose from a normal centered at current value

proposal <- rnorm(1, mean = mu, sd = 0.5)

# Log posterior (up to constant): -n/2 _ (ybar - mu)^2 - mu^2 / (2 _ 100)

log*post <- function(m) {
-n/2 * (ybar - m)^2 - m^2 / (2 _ 100)
}

# Accept ratio

log_R <- log_post(proposal) - log_post(mu)

if (log(runif(1)) < log_R) {
mu <- proposal
}
samples[t] <- mu
}

cat("Metropolis-Hastings results:\n")
cat(" Posterior mean:", round(mean(samples), 3), "\n")
cat(" Posterior SD:", round(sd(samples), 3), "\n")
cat(" True posterior mean:", round(ybar, 3), "\n")

# Trace plot

plot(samples, type = "l", col = "steelblue",
main = "MCMC trace plot", xlab = "Iteration", ylab = "mu")

# Histogram

hist(samples, breaks = 50, col = "steelblue", border = "white",
main = "MCMC posterior samples", xlab = "mu", freq = FALSE)
curve(dnorm(x, mean = ybar, sd = 1/sqrt(n)), add = TRUE, lwd = 2, col = "firebrick")
legend("topright", legend = c("MCMC", "Exact"), col = c("steelblue", "firebrick"), lwd = 2)

Click Run to execute this code in your browser with WebR

8.2.3 Muestreo de Gibbs

El muestreo de Gibbs es un caso especial de Metropolis-Hastings donde cada parámetro se actualiza condicionalmente a todos los demás. Es eficiente cuando las distribuciones condicionales son fáciles de muestrear.

Rethinking

MCMC es una herramienta, no una caja negra. La belleza de MCMC es que funciona para cualquier distribución a posteriori. Pero este poder viene con responsabilidad: debe verificar que la cadena ha convergido. MCMC puede quedarse atascado, explorar lentamente, o perder regiones de la distribución a posteriori. Los diagnósticos (Sección 8.4) ayudan a detectar estos problemas. La clave es entender la mecánica — no tratar MCMC como una caja mágica que produce la respuesta.

Resumen de la Sección 8.2

MCMC muestrea de la distribución a posteriori sin conocer la constante de normalización.
Funciona para cualquier forma de distribución a posteriori, a diferencia de la aproximación cuadrática.
Metropolis-Hastings generaliza el algoritmo de Metropolis con propuestas arbitrarias.
El muestreo de Gibbs actualiza cada parámetro condicionalmente a los demás.
MCMC es una herramienta, no una caja negra — siempre verifique la convergencia.

El buen rey Markov y su reino insular

HMC fácil: map2stan