Algoritmo Metropolis–Hastings

Introducción

El algoritmo Metropolis–Hastings (MH) es un método de simulación Monte Carlo basado en cadenas de Markov (MCMC) utilizado para generar muestras de una distribución de probabilidad cuando su forma analítica es difícil o imposible de obtener.

En este documento se ilustra el funcionamiento del algoritmo aplicándolo a una distribución Normal(3,1), mostrando cómo las muestras simuladas se ajustan a la densidad teórica esperada.

---

Marco teórico

En inferencia bayesiana, el algoritmo Metropolis–Hastings permite generar una cadena de valores que convergen a la distribución posterior objetivo \(\pi(\theta)\).
La probabilidad de aceptación se define como:

\[ \alpha = \min\left(1, \frac{\pi(\theta^*)}{\pi(\theta^{(t-1)})}\right) \]

Cuando la propuesta es simétrica (por ejemplo, una normal centrada en el valor actual), esta forma es suficiente para obtener la distribución deseada.

---

Implementación en R

A continuación se muestra la implementación paso a paso del algoritmo Metropolis–Hastings en R, con una distribución objetivo \(N(3,1)\).

set.seed(123)      
n <- 10000         # Número de iteraciones
theta <- numeric(n)
theta[1] <- 0      # Valor inicial
sigma_prop <- 1    # Desviación estándar de la propuesta

# Función objetivo (densidad proporcional)
target <- function(x) {
  exp(-0.5 * (x - 3)^2)
}

# Bucle principal del algoritmo
for (t in 2:n) {
  proposal <- rnorm(1, mean = theta[t-1], sd = sigma_prop)
  alpha <- min(1, target(proposal) / target(theta[t-1]))
  
  if (runif(1) < alpha) {
    theta[t] <- proposal
  } else {
    theta[t] <- theta[t-1]
  }
}

# Muestras obtenidas
summary(theta)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.3541  2.3179  3.0093  2.9954  3.6761  6.8948 

Resultados

A continuación se presentan las dos gráficas principales que evidencian el funcionamiento del algoritmo.

Histograma y densidad teórica

x <- seq(0, 6, length.out = 200)
dens_teorica <- dnorm(x, mean = 3, sd = 1)

hist(theta, breaks = 40, freq = FALSE, col = "lightblue",
main = "Muestras del Algoritmo Metropolis–Hastings",
xlab = expression(theta), ylab = "Densidad")
lines(x, dens_teorica, col = "red", lwd = 2)
legend("topright", legend = c("Muestras MH", "Normal(3,1)"),
col = c("lightblue", "red"), lwd = c(10, 2), bty = "n")

Interpretación

Las muestras simuladas siguen la forma teórica de la distribución Normal(3,1), lo que indica que el algoritmo explora correctamente la región de mayor probabilidad y converge a la distribución objetivo.

Trace plot de las primeras 500 iteraciones

plot(theta[1:500], type = "l", lwd = 0.8,
main = "Trace plot del algoritmo Metropolis–Hastings",
xlab = "Iteración", ylab = expression(theta))

Interpretación

El gráfico muestra oscilaciones alrededor del valor esperado (\(\theta = 3\)), sin tendencias crecientes ni decrecientes. Esto significa que la cadena alcanzó el equilibrio y que el muestreo es estable.

Trace plot completo

plot(theta, type = "l", lwd = 0.8,
main = "Trace plot del algoritmo Metropolis–Hastings",
xlab = "Iteración", ylab = expression(theta))