8.6 Práctica

Los problemas de práctica del Capítulo 8 le piden ejecutar MCMC en modelos simples, diagnosticar la convergencia y comparar HMC con aproximación por rejilla.

Practice

Problemas fáciles

8E1. Dibuje un DAG para el algoritmo de Metropolis. ¿Cuáles son los nodos y cuáles las flechas?

8E2. ¿Por qué el algoritmo de Metropolis solo necesita el cociente de probabilidades, no los valores absolutos?

8E3. ¿Cuál es la diferencia entre Metropolis, Metropolis-Hastings y muestreo de Gibbs?

8E4. ¿Qué mide R-hat? ¿Qué valor indica convergencia?

Problemas medios

8M1. Implemente el algoritmo de Metropolis para una distribución a posteriori simple y compare las muestras con la distribución a posteriori exacta.

# 8M1: Metropolis algorithm for a simple posterior
# Model: y ~ Normal(mu, 1), prior: mu ~ Normal(0, 10)
# Exact posterior: mu ~ Normal(ybar, 1/n)

set.seed(42)
n <- 50
y <- rnorm(n, mean = 5, sd = 1)
ybar <- mean(y)

# Metropolis

n_iter <- 10000
mu <- 0
samples <- numeric(n_iter)

for (t in 1:n_iter) {
proposal <- mu + rnorm(1, sd = 0.3)
log_post <- function(m) -n/2 * (ybar - m)^2 - m^2 / 200
if (log(runif(1)) < log_post(proposal) - log_post(mu)) mu <- proposal
samples[t] <- mu
}

cat("Metropolis:\n")
cat(" Mean:", round(mean(samples), 3), "\n")
cat(" SD:", round(sd(samples), 3), "\n")

cat("\nExact posterior:\n")
cat(" Mean:", round(ybar, 3), "\n")
cat(" SD:", round(1/sqrt(n), 3), "\n")

# Compare histograms

hist(samples, breaks = 50, col = adjustcolor("steelblue", 0.5),
freq = FALSE, main = "8M1: Metropolis vs exact posterior",
xlab = "mu")
curve(dnorm(x, mean = ybar, sd = 1/sqrt(n)), add = TRUE, lwd = 3, col = "firebrick")
legend("topright", legend = c("Metropolis", "Exact"), col = c("steelblue", "firebrick"), lwd = 2)

Click Run to execute this code in your browser with WebR

8M2. Ejecute el algoritmo de Metropolis con diferentes desviaciones estándar de propuesta (0.1, 0.5, 1.0). ¿Cómo afecta la propuesta a la cadena?

# 8M2: Effect of proposal SD on Metropolis
set.seed(42)
n <- 50
y <- rnorm(n, mean = 5, sd = 1)
ybar <- mean(y)

log_post <- function(m) -n/2 * (ybar - m)^2 - m^2 / 200

run_metropolis <- function(prop_sd, n_iter = 3000) {
mu <- 0
samples <- numeric(n_iter)
accepts <- 0
for (t in 1:n_iter) {
proposal <- mu + rnorm(1, sd = prop_sd)
if (log(runif(1)) < log_post(proposal) - log_post(mu)) {
mu <- proposal
accepts <- accepts + 1
}
samples[t] <- mu
}
return(list(samples = samples, accept_rate = accepts / n_iter))
}

par(mfrow = c(3, 1))
for (sd in c(0.1, 0.5, 1.0)) {
result <- run_metropolis(sd)
plot(result$samples, type = "l", col = "steelblue",
     main = sprintf("Proposal SD = %.1f (accept rate = %.2f)", sd, result$accept_rate),
xlab = "Iteration", ylab = "mu")
}
par(mfrow = c(1, 1))

cat("Small proposal: high accept rate, slow exploration\n")
cat("Large proposal: low accept rate, fast but jumpy\n")
cat("Medium proposal: balanced accept rate and exploration\n")

Click Run to execute this code in your browser with WebR

Problemas difíciles

8H1. Implemente un muestreador HMC simple para una distribución a posteriori gaussiana bidimensional. Verifique la convergencia con R-hat.

# 8H1: Simple HMC for a 2D Gaussian
set.seed(42)

# Target: 2D standard normal

log_post <- function(theta) -0.5 * sum(theta^2)
grad_log_post <- function(theta) -theta

# HMC sampler

hmc_2d <- function(n_iter, step_size = 0.2, n_leapfrog = 15) {
theta <- c(0, 0) # starting value
samples <- matrix(0, nrow = n_iter, ncol = 2)

for (t in 1:n_iter) { # Random momentum
p <- rnorm(2)

  # Leapfrog
  theta_new <- theta
  p_new <- p
  for (i in 1:n_leapfrog) {
    p_new <- p_new + 0.5 * step_size * grad_log_post(theta_new)
    theta_new <- theta_new + step_size * p_new
    p_new <- p_new + 0.5 * step_size * grad_log_post(theta_new)
  }

  # Metropolis accept
  current_H <- -log_post(theta) + 0.5 * sum(p^2)
  proposed_H <- -log_post(theta_new) + 0.5 * sum(p_new^2)
  if (log(runif(1)) < current_H - proposed_H) theta <- theta_new
  samples[t, ] <- theta

}
return(samples)
}

# Run 4 chains

chains <- list()
for (c in 1:4) {
chains[[c]] <- hmc_2d(2000)
}

# Compute R-hat for each parameter

compute*rhat <- function(chain_list) {
m <- length(chain_list)
n <- nrow(chain_list[[1]])
mat <- sapply(chain_list, function(x) x[, 1])
W <- mean(apply(mat, 2, var))
B <- n * var(colMeans(mat))
var*hat <- (n - 1) / n * W + B / n
return(sqrt(var_hat / W))
}

cat("HMC 2D Gaussian:\n")
cat(" Parameter 1: mean =", round(mean(sapply(chains, function(x) mean(x[,1]))), 3),
" R-hat =", round(compute_rhat(chains), 4), "\n")

# Plot

plot(chains[[1]][,1], chains[[1]][,2], type = "l", col = "steelblue",
main = "8H1: HMC 2D Gaussian", xlab = "theta1", ylab = "theta2")
for (c in 2:4) {
lines(chains[[c]][,1], chains[[c]][,2], col = "steelblue")
}
points(0, 0, pch = 16, col = "firebrick", cex = 2)
legend("topright", legend = c("HMC samples", "True mean"),
col = c("steelblue", "firebrick"), pch = c(NA, 16), lwd = c(1, NA))

Click Run to execute this code in your browser with WebR

8H2. Compare el tamaño efectivo de muestra del Metropolis de caminata aleatoria y HMC para la misma distribución objetivo.

8H3. Simule datos de un modelo con una distribución a posteriori no normal (ej., una mezcla). Ajuste con quap y con HMC. ¿Cómo difieren los resultados?

Resumen

Entropía máxima