Metropolis Method

Metropolis Method

El algoritmo de :contentReferenceoaicite:0 es un método de simulación basado en cadenas de Markov que permite generar muestras de una distribución objetivo \(\pi(x)\) cuando no es posible muestrear directamente de ella. Este algoritmo es fundamental dentro de los métodos de Monte Carlo vía cadenas de Markov (MCMC) y tiene amplias aplicaciones en inferencia bayesiana.

Idea general

El objetivo es construir una cadena de Markov \(\{X_t\}\) cuya distribución estacionaria sea \(\pi(x)\). Para ello, se utiliza una distribución propuesta \(q(x^*|x)\) que sugiere nuevos estados candidatos.

En cada iteración, el algoritmo decide si acepta o rechaza el nuevo valor basándose en una probabilidad que garantiza la convergencia hacia la distribución objetivo.

Regla de aceptación

La probabilidad de aceptación está dada por:

\[ \alpha(x_t, x^*) = \min\left(1, \frac{\pi(x^*) q(x_t|x^*)}{\pi(x_t) q(x^*|x_t)} \right) \]

Esta expresión asegura que la cadena satisface la propiedad de equilibrio detallado, lo que implica que \(\pi(x)\) es distribución estacionaria.

Algoritmo

1. Inicializar \(x_0\)
2. Para \(t = 1,2,\dots,n\):
   • Generar \(x^* \sim q(x^*|x_t)\)
   • Calcular \(\alpha(x_t, x^*)\)
   • Generar \(u \sim U(0,1)\)
   • Si \(u < \alpha\), aceptar: \(x_{t+1} = x^*\)
   • En caso contrario, rechazar: \(x_{t+1} = x_t\)

Caso especial: algoritmo de Metropolis

Cuando la distribución propuesta es simétrica, es decir,

\[ q(x^*|x) = q(x|x^*) \]

la probabilidad de aceptación se simplifica a:

\[ \alpha(x_t, x^*) = \min\left(1, \frac{\pi(x^*)}{\pi(x_t)} \right) \]

Propiedades

• No requiere conocer la constante de normalización de \(\pi(x)\)
• Genera muestras dependientes pero convergentes
• Es válido en espacios de alta dimensión
• Depende críticamente de la elección de la propuesta \(q\)

Ejemplo en R

set.seed(123)

# Distribución objetivo (no normalizada)
pi_target <- function(x) {
  exp(-x^2 / 2)
}

n <- 10000
x <- numeric(n)
x[1] <- 0
sigma <- 1

for (t in 2:n) {
  x_star <- rnorm(1, mean = x[t-1], sd = sigma)
  alpha <- min(1, pi_target(x_star) / pi_target(x[t-1]))
  
  if (runif(1) < alpha) {
    x[t] <- x_star
  } else {
    x[t] <- x[t-1]
  }
}

hist(x, breaks = 50, probability = TRUE, col = "lightblue",
     main = "Simulación Metropolis-Hastings")
curve(dnorm(x), add = TRUE, col = "red", lwd = 2)

set.seed(123)

# Distribución objetivo (no normalizada)
pi_target <- function(x) {
  exp(-x^2 / 2)
}

# Parámetros
n <- 10000
x <- numeric(n)
x[1] <- 0
sigma <- 1  # varianza de propuesta

for (t in 2:n) {
  # Propuesta normal centrada en el estado actual
  x_star <- rnorm(1, mean = x[t-1], sd = sigma)
  
  # Probabilidad de aceptación
  alpha <- min(1, pi_target(x_star) / pi_target(x[t-1]))
  
  # Aceptar o rechazar
  if (runif(1) < alpha) {
    x[t] <- x_star
  } else {
    x[t] <- x[t-1]
  }
}

# Visualización
hist(x, breaks = 50, probability = TRUE, col = "lightblue")
curve(dnorm(x), add = TRUE, col = "red", lwd = 2)

set.seed(123)

# Densidad objetivo (no normalizada)
pi_target <- function(x) {
  ifelse(x >= 0 & x <= 3, x^2 + x^3 + x^8, 0)
}

# Parámetros
n <- 20000
x <- numeric(n)
x[1] <- 1  # valor inicial
sigma <- 0.5

for (t in 2:n) {
  # Propuesta
  x_star <- rnorm(1, mean = x[t-1], sd = sigma)
  
  # Si cae fuera del soporte, se rechaza automáticamente
  if (x_star < 0 || x_star > 3) {
    x[t] <- x[t-1]
  } else {
    # Probabilidad de aceptación
    alpha <- min(1, pi_target(x_star) / pi_target(x[t-1]))
    
    if (runif(1) < alpha) {
      x[t] <- x_star
    } else {
      x[t] <- x[t-1]
    }
  }
}

# Burn-in
burn <- 2000
x_sample <- x[(burn+1):n]

# Histograma
hist(x_sample, probability = TRUE, breaks = 50,
     col = "lightblue", main = "Muestra Metropolis-Hastings")

# Curva teórica (normalizada numéricamente)
f <- function(x) x^2 + x^3 + x^8
Z <- integrate(f, lower = 0, upper = 3)$value

curve(f(x)/Z, from = 0, to = 3, col = "red", lwd = 2, add = TRUE)