Modelo

A continuación se describe la estructura básica de un modelo de Regresión Poisson como un modelo lineal generalizado (GLM, Generalized Linear Model).

Distribución muestral:

\[ y_i \mid \theta_i \stackrel{\text{iid}}{\sim} \textsf{Poisson}(\theta_i), \qquad i = 1, \ldots, n, \] donde la media \(\theta_i\) está relacionada con una combinación lineal de predictores a través de una transformación logarítmica: \[ \log(\theta_i) = \boldsymbol{\beta}^{\textsf{T}} \boldsymbol{x}_i = \sum_{j=1}^p \beta_j x_{i,j}, \] donde \(\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)\) y \(\boldsymbol{x}i = (x{i1}, \ldots, x_{ip})\) denotan, respectivamente, el vector de coeficientes desconocidos y el vector de predictores observados correspondientes a la \(i\)-ésima observación.

La función de enlace (link function) utilizada es la logarítmica, \(\log\), la cual asegura que \(\theta_i > 0\), en concordancia con las propiedades de la distribución Poisson.

Distribución previa:

No existen distribuciones previas conjugadas para los modelos lineales generalizados, excepto en el caso de la regresión normal. Por ello, se emplea una distribución previa normal como aproximación: \[ \boldsymbol{\beta} \sim \textsf{N}_p(\boldsymbol{\beta}_0, \Sigma_0)\,, \] donde \(\boldsymbol{\beta}_0\) es la media previa de los coeficientes y \(\Sigma_0\) esl la matriz de covarianza previa.

Los parámetros del modelo son \((\beta_1, \ldots, \beta_p)\).

Los hiperparámetros del modelo son \((\boldsymbol{\beta}_0, \Sigma_0)\).

Distibución conjunta posterior

(Ejercicio.) La distribución posterior del modelo es de la forma: \[ p(\boldsymbol{\beta}\mid\boldsymbol{y}) \propto \prod_{i=1}^n e^{-\theta_i}\,\theta_i^{y_i} \times \exp\left\{ -\tfrac{1}{2} \left[ \boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta} - 2\boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta}_0 \right] \right\}\,, \] donde \(\log(\theta_i) = \boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i\).

Observe que \(p(\boldsymbol{\beta}\mid\boldsymbol{y})\) no corresponde a ninguna familia paramétrica de distribuciones conocida.

(Ejercicio.) En escala logarítmica, se tiene que \[ \log p(\boldsymbol{\beta}\mid\boldsymbol{y}) = \sum_{i=1}^n(y_i\log(\theta_i) - \theta_i) -\tfrac{1}{2} \left[ \boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta} - 2\boldsymbol{\beta}^{\textsf{T}}\Sigma_0^{-1}\boldsymbol{\beta}_0 \right] + \text{C} \] donde \(\text{C}\) es una constante (términos que no dependen de \(\boldsymbol{\beta}\)).

(Ejercicio.) El gradiente correspondiente está dado por \[ \frac{\partial}{\partial \boldsymbol{\beta}}\log p(\boldsymbol{\beta} \mid \boldsymbol{y}) = \sum_{i=1}^n \left( y_i - \exp{\left(\boldsymbol{\beta}^{\textsf{T}}\boldsymbol{x}_i\right)} \right)\boldsymbol{x}_i\,. \]

Algorithmo de Metropolis-Hastings

El algoritmo de Metropolis-Hastings es un método general para aproximar distribuciones probabilísticas, incluso cuando no se dispone de una solución analítica exacta.

El principal desafío aparece cuando \(p(\boldsymbol{\theta} \mid \boldsymbol{y})\) no corresponde a una distribución estándar conocida, lo que impide generar muestras directamente de la misma.

Algoritmo de Metropolis

El algoritmo de Metropolis es un caso particular del algoritmo de Metropolis-Hastings.

Este algoritmo produce una cadena de Markov que converge a la distribución posterior objetivo \(p(\theta \mid \boldsymbol{y})\), bajo el supuesto de que la distribución de propuestas es simétrica.

Algoritmo

Dado el estado actual \(\theta^{(b)}\) del parámetro \(\theta\), el siguiente estado \(\theta^{(b+1)}\) se genera de la siguiente manera:

Generación de la propuesta:
Simular un valor propuesto \(\theta^* \sim J(\theta \mid \theta^{(b)})\), donde \(J\) es una distribución simétrica tal que \(J(\theta_1 \mid \theta_2) = J(\theta_2 \mid \theta_1)\).
Esta distribución, conocida como distribución de propuestas (proposal distribution), suele especificarse como:
\[ J(\theta \mid \theta^{(b)}) = \textsf{N}(\theta \mid \theta^{(b)}, \delta^2) \quad \text{o} \quad J(\theta \mid \theta^{(b)}) = \textsf{U}(\theta \mid \theta^{(b)} - \delta, \theta^{(b)} + \delta), \] donde \(\delta\) es el parámetro de ajuste, elegido para garantizar un rendimiento adecuado del algoritmo.
Cálculo de la probabilidad de aceptación:
Calcular la probabilidad de aceptación \(r\), definida como:
\[ r = \frac{p(\theta^* \mid \boldsymbol{y})}{p(\theta^{(b)} \mid \boldsymbol{y})}. \]
Para mejorar la estabilidad numérica, este cálculo suele realizarse en la escala logarítmica:
\[ r = \exp\left[ \log p(\theta^* \mid \boldsymbol{y}) - \log p(\theta^{(b)} \mid \boldsymbol{y}) \right]. \]
Actualización del estado:
Actualizar el estado \(\theta^{(b+1)}\) según la probabilidad de aceptación:
\[ \theta^{(b+1)} = \begin{cases} \theta^*, & \text{con probabilidad } \min(1, r), \\ \theta^{(b)}, & \text{con probabilidad } 1 - \min(1, r). \end{cases} \]

Consideraciones

Este algoritmo produce una secuencia de valores \(\theta^{(1)}, \ldots, \theta^{(B)}\), donde \(\theta^{(b+1)}\) depende únicamente del estado anterior \(\theta^{(b)}\). Esto hace que la secuencia sea una cadena de Markov, lo que garantiza su convergencia a la distribución objetivo bajo condiciones apropiadas.
Es común calibrar \(\delta\) para lograr una tasa de aceptación entre el 30% y el 50% (0.234 para cadenas de alta dimensión o 0.44 para dimensión 1), lo que suele optimizar el rendimiento del algoritmo. Para facilitar esta calibración, pueden emplearse algoritmos adaptativos que ajusten automáticamente el valor de \(\delta\) durante las iteraciones iniciales.
El parámetro de ajuste \(\delta\) controla el tamaño de los pasos propuestos y tiene un impacto directo en la eficiencia del algoritmo. Valores de \(\delta\) demasiado pequeños generan una exploración lenta del espacio de parámetros, mientras que valores grandes pueden conducir a bajas tasas de aceptación.
La correlación de la cadena entre estados consecutivos depende del valor de \(\delta\) en la distribución de propuestas. Una correlación alta indica una exploración menos eficiente del espacio de parámetros.

Ejemplo: Modelo Normal

Algoritmo de Metropolis en un modelo Normal con varianza conocida.

Los datos \(y_i\) se suponen independientes e idénticamente distribuidos (i.i.d.) según una distribución Normal con media \(\theta = 10\) y varianza \(1\).

El modelo es de la forma: \[ y_i \mid \theta \overset{\text{i.i.d.}}{\sim} \textsf{N}(\theta, 1), \qquad i = 1, \ldots, n. \]

Se especifica una distribución previa conjugada para el parámetro \(\theta\):
\[ \theta \sim \textsf{N}(\mu, \tau^2). \]

# Reproducibilidad
set.seed(123)

# Simular datos
n  <- 5
s2 <- 1
y  <- rnorm(n = n, mean = 10, sd = 1)

# Previa
t2 <- 10
mu <- 5

# Algoritmo de Metropolis
theta <- 0                              # valor inicial
delta <- 1.75                           # parámetro de ajuste
acr   <- 0                              # tasa de aceptación
B     <- 10000                          # número de muestras
THETA <- matrix(NA, nrow = B, ncol = 1)  # almacenamiento

# cadena
for (b in 1:B) {
  # 1. propuesta
  theta_star <- rnorm(n = 1, mean = theta, sd = sqrt(delta))

  # 2. probabilidad de aceptación (en escala log)
  log_post_star <- sum(dnorm(x = y, mean = theta_star, sd = sqrt(s2), log = TRUE)) +
                   dnorm(x = theta_star, mean = mu, sd = sqrt(t2), log = TRUE)

  log_post_curr <- sum(dnorm(x = y, mean = theta, sd = sqrt(s2), log = TRUE)) +
                   dnorm(x = theta, mean = mu, sd = sqrt(t2), log = TRUE)

  log_r <- log_post_star - log_post_curr

  # 3. actualizar
  if (runif(1) < exp(log_r)) {
    theta <- theta_star
    acr <- acr + 1/B
  }

  # 4. almacenar
  THETA[b] <- theta
}

# Tasa de aceptación
round(100 * acr, 1)

## [1] 38

# Gráfico
par(mfrow = c(1, 2), mar = c(3, 3, 1, 1), mgp = c(1.75, 0.75, 0))

# Cadena
plot(x = seq_along(THETA), y = THETA, type = "p", pch = 20, cex = 0.3,
  xlab = "Iteración", ylab = expression(theta))

# Histograma (sin periodo de calentamiento)
hist(x = THETA[-(1:10)], freq = FALSE, col = "gray", border = "gray",
  main = "", xlab = expression(theta), ylab = "Densidad")

# Posterior analítica
th   <- seq(min(THETA), max(THETA), length.out = 1000)
mu_n <- (mu / t2 + n * mean(y) / s2) / (1 / t2 + n / s2)
t2_n <- 1 / (1 / t2 + n / s2)

lines(x = th, y = dnorm(th, mean = mu_n, sd = sqrt(t2_n)),
      col = 2, lty = 1, lwd = 2)

# Reproducibilidad
set.seed(123)

ACR <- NULL  # tasa de aceptaciones
ACF <- NULL  # autocorrelaciones
SAM <- NULL  # muestras

for (delta2 in 2^c(-5, -1, 1, 5, 7)) {
  # parámetros iniciales
  theta <- 0                               # valor inicial
  acr   <- 0                               # tasa de aceptación
  B     <- 10000                           # número de muestras
  THETA <- matrix(NA, nrow = B, ncol = 1)  # almacenamiento

  # cadena
  for (b in 1:B) {
    # 1. propuesta
    theta_star <- rnorm(n = 1, mean = theta, sd = sqrt(delta2))

    # 2. probabilidad de aceptación (escala log)
    log_post_star <- sum(dnorm(x = y, mean = theta_star, sd = sqrt(s2), log = TRUE)) +
                   dnorm(x = theta_star, mean = mu, sd = sqrt(t2), log = TRUE)

    log_post_curr <- sum(dnorm(x = y, mean = theta, sd = sqrt(s2), log = TRUE)) +
                   dnorm(x = theta, mean = mu, sd = sqrt(t2), log = TRUE)

    log_r <- log_post_star - log_post_curr

    # 3. actualizar
    if (runif(1) < exp(log_r)) {
      theta <- theta_star
      acr <- acr + 1/B
    }

    # 4. almacenar
    THETA[b] <- theta
  }

  # almacenar valores para todos los valores de delta2
  ACR <- c(ACR, acr)
  ACF <- c(ACF, acf(THETA, plot = FALSE)$acf[2])
  SAM <- cbind(SAM, THETA)
}

# Parámetros de ajuste
round(2^c(-5, -1, 1, 5, 7), 3)

## [1]   0.031   0.500   2.000  32.000 128.000

# Tasas de aceptación
round(100 * ACR, 1)

## [1] 86.4 56.3 35.4  8.8  4.8

# Autocorrelaciones
round(ACF, 3)

## [1] 0.975 0.845 0.743 0.891 0.915

# Tamaños efectivos de muestra
round(coda::effectiveSize(SAM), 0)

## var1 var2 var3 var4 var5 
##  120  347 1027  611  399

Algoritmo de Metropolis adaptativo

Emplear la varianza posterior en la distribución de propuestas suele ser la opción más eficiente.

Aunque la varianza posterior exacta no se conoce antes de ejecutar el algoritmo de Metropolis, una aproximación suele ser suficiente para lograr un buen desempeño.

Si la tasa de aceptación obtenida es demasiado baja o demasiado alta, puede ajustarse la variabilidad de la propuesta de forma adaptativa con el fin de optimizar el rendimiento del algoritmo.

Una forma común de adaptar el parámetro de ajuste (tuning parameter) \(\delta\) en el algoritmo de Metropolis es usar una regla de actualización estocástica basada en la tasa de aceptación observada: \[ \log \delta_{t+1} \;=\; \log \delta_t \;+\; \gamma_t \,(\alpha_t - \alpha^\ast) \qquad\Longleftrightarrow\qquad \delta_{t+1} \;=\; \delta_t \; \exp\!\big(\,\gamma_t \,(\alpha_t - \alpha^\ast)\,\big), \] donde:

\(\delta_t\): parámetro de ajuste en la iteración \(t\).
\(\alpha_t\): tasa de aceptación observada en la iteración \(t\).
\(\alpha^\ast\): tasa de aceptación objetivo (0.234 para cadenas de alta dimensión o 0.44 para dimensión 1).
\(\gamma_t\): secuencia de factores de paso (step sizes) tales que \(\sum_t \gamma_t = \infty\) y \(\sum_t \gamma_t^2 < \infty\) (e.g., \(\gamma_t = t^{-0.5}\)).

¡Cuidado! La adaptación se realiza solo durante el periodo de calentamiento (burn-in).

Algoritmo Metropolis-Hastings

El algoritmo Metropolis-Hastings es una generalización del algoritmo de Metropolis, que permite utilizar una distribución de propuestas arbitraria, lo que lo hace más flexible para aplicaciones donde estas variantes más simples son inadecuadas.

Algoritmo

Dado el estado actual \(\theta^{(b)}\) del parámetro \(\theta\), el siguiente estado \(\theta^{(b+1)}\) se obtiene mediante el siguiente procedimiento:

Generación de la propuesta:
Simular un valor propuesto \(\theta^* \sim J(\theta \mid \theta^{(b)})\), donde \(J\) es una distribución de propuestas arbitraria.
Cálculo de la tasa de aceptación:
Evaluar la tasa de aceptación \(r\) como
\[ r = \frac{p(\theta^* \mid y)}{p(\theta^{(b)} \mid y)} \cdot \frac{J(\theta^{(b)} \mid \theta^*)}{J(\theta^* \mid \theta^{(b)})}. \]
Este ajuste incorpora el cociente de las probabilidades de la distribución de propuestas para corregir posibles asimetrías en \(J\).
Actualización del estado:
Definir \(\theta^{(b+1)}\) de acuerdo con la probabilidad de aceptación:
\[ \theta^{(b+1)} = \begin{cases} \theta^*, & \text{con probabilidad } \min(1, r), \\ \theta^{(b)}, & \text{con probabilidad } 1 - \min(1, r). \end{cases} \]

Consideraciones

El algoritmo de Metropolis es un caso particular del método de Metropolis–Hastings en el que la distribución de propuestas \(J(\theta \mid \theta^{(b)})\) es simétrica, es decir, \(J(\theta_1 \mid \theta_2) = J(\theta_2 \mid \theta_1)\).
La simetría de \(J(\theta \mid \theta^{(b)})\) simplifica el cálculo de la probabilidad de aceptación, pues elimina la necesidad de incluir el cociente de probabilidades de la distribución de propuestas, reduciendo así la complejidad computacional.

Ejemplo: Regresión Poisson

Investigar las actividades de reproducción de gorriones en función de la edad.

El estudio analiza la relación entre la edad de los gorriones y el número de crías que producen, con el objetivo de caracterizar cómo la edad influye en su éxito reproductivo.

Referencias:

Arcese, P., Smith, J. N. M., Hochachka, W. M., Rogers, C. M., & Ludwig, D. (1992). Stability, regulation, and the determination of abundance in an insular song sparrow population. Ecology, 73(3), 805–822.
Hoff, P. D. (2009). A first course in Bayesian statistical methods. New York: Springer.

Tratamiento de datos

Las variables del modelo son:

\(y_i\): número de crías producidas por el individuo \(i\), para \(i = 1, \ldots, n\).
\(x_{i,j}\): valor del predictor \(j\) observado en el individuo \(i\), para \(i = 1, \ldots, n\) y \(j = 1, \ldots, k\).

# Actividades de reproducción de gorriones (hembras) en función de la edad
# "age"     : edad.
# "fledged" : número de crías.
dat <- structure(
  c(
    3, 1, 1, 2, 0, 0, 6, 3, 4, 2, 1, 6, 2, 3, 3, 4, 7, 2, 2, 1, 
    1, 3, 5, 5, 0, 2, 1, 2, 6, 6, 2, 2, 0, 2, 4, 1, 2, 5, 1, 2, 
    1, 0, 0, 2, 4, 2, 2, 2, 2, 0, 3, 2, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 
    2, 2, 2, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 4, 4, 
    4, 4, 5, 5, 5, 5, 3, 3, 3, 3, 3, 3, 3, 6, 1, 1, 9, 9, 1, 1, 
    1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 25, 25, 16, 16, 
    16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 25, 16, 16, 16, 16, 
    25, 25, 25, 25, 9, 9, 9, 9, 9, 9, 9, 36, 1, 1
  ),
  .Dim = c(52L, 4L),
  .Dimnames = list(NULL, c("fledged", "intercept", "age", "age2"))
)
dat <- as.data.frame(dat)

y    <- dat$fledged  # y  = variable dependiente (respuesta)
age  <- dat$age      # x1 = variable independiente 1
age2 <- age^2        # x2 = variable independiente 2

Análisis exploratorio

# Tendencia
summary(y)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   2.000   2.404   3.000   7.000

# Variabilidad
round(var(y), 3)

## [1] 3.265

# Visualización
par(mar = c(3, 3, 1, 1), mgp = c(1.75, 0.75, 0))

# Boxplots
plot(y ~ as.factor(age), boxwex = 0.4, range = 0,
  xlab = "Edad (años)", ylab = "No. Crias", col = 0, border = "darkgray")

# Añadir observaciones con jitter
points(jitter(as.numeric(as.factor(age)), amount = 0.15), y,
  pch = 20, col = adjustcolor(4, 0.5), cex = 1.2)

# GLM (frecuentista)
summary(glm(y ~ age + age2, family="poisson"))

## 
## Call:
## glm(formula = y ~ age + age2, family = "poisson")
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.27662    0.44219   0.626   0.5316  
## age          0.68174    0.33850   2.014   0.0440 *
## age2        -0.13451    0.05786  -2.325   0.0201 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 76.081  on 51  degrees of freedom
## Residual deviance: 67.837  on 49  degrees of freedom
## AIC: 198.78
## 
## Number of Fisher Scoring iterations: 5

Ajuste del modelo

# Variable respuesta
y <- dat$fledged

# Matriz de diseño (incluye intercepto, edad y edad^2)
X <- cbind(1, age, age^2)

# Tamaño de la muestra y número de predictores
(n <- nrow(X))

## [1] 52

(p <- ncol(X))

## [1] 3

# Hiperparámetros
beta0  <- rep(0, p)   # beta0 = 0
Sigma0 <- rep(10, p)  # Sigma0 = 10*I

# Parámetro de ajuste
Delta <- var(log(y + 1)) * solve(t(X)%*%X)

# Reproducibilidad
set.seed(123)

# Algoritmo de Metropolis
beta <- rep(0, p)                       # valor inicial
acr  <- 0                               # tasa de aceptación
B    <- 100000                          # número de muestras
BETA <- matrix(NA, nrow = B, ncol = p)  # almacenamiento

# Progreso
ncat <- floor(B / 10)

# Cadena
for (b in 1:B) {
  # 1. Propuesta
  beta_p <- c(mvtnorm::rmvnorm(n = 1, mean = beta, sigma = Delta))

  # 2. Tasa de aceptación (escala log)
  eta_p <- as.numeric(X %*% beta_p)    # predictor lineal propuesto
  eta   <- as.numeric(X %*% beta)      # predictor lineal actual

  log_post_star <- sum(dpois(x = y, lambda = exp(eta_p), log = TRUE)) + 
       sum(dnorm(x = beta_p, mean = beta0, sd = Sigma0, log = TRUE))
       
  log_post_curr <- sum(dpois(x = y, lambda = exp(eta), log = TRUE)) + 
       sum(dnorm(x = beta, mean = beta0, sd = Sigma0, log = TRUE))

  log_r <- log_post_star - log_post_curr

  # 3) Actualizar
  if (runif(1) < exp(log_r)) {
    beta <- beta_p
    acr  <- acr + 1/B
  }

  # 4. Almacenar
  BETA[b, ] <- beta

  # 5. Progreso
  if (b %% ncat == 0) {
    cat(round(100 * b / B, 1), "% completado ...\n", sep = "")
  }
}

## 10% completado ...
## 20% completado ...
## 30% completado ...
## 40% completado ...
## 50% completado ...
## 60% completado ...
## 70% completado ...
## 80% completado ...
## 90% completado ...
## 100% completado ...

# Tasa de aceptación
round(100 * acr, 1)

## [1] 53

# Tamaños efectivos de muestra
round(apply(X = BETA, MARGIN = 2, FUN = coda::effectiveSize))

## [1] 6379 6008 5524

# Intervalos de credibilidad al 95%
(round(quantile(BETA[, 2], probs = c(0.025, 0.975)), 3))

##  2.5% 97.5% 
## 0.073 1.408

(round(quantile(BETA[, 3], probs = c(0.025, 0.975)), 3))

##   2.5%  97.5% 
## -0.259 -0.032

# Probabilidades posteriores de ser > 0
(round(mean(BETA[, 2] > 0), 3))

## [1] 0.986

(round(mean(BETA[, 3] > 0), 3))

## [1] 0.005

La inferencia sobre \(\beta_2\) indica que el número de crías aumenta con la edad en las primeras etapas de vida, mientras que la estimación de \(\beta_3\) revela un efecto cuadrático negativo, lo que implica que el éxito reproductivo declina en edades avanzadas.

Algoritmo de Monte Carlo Hamiltoniano

El algoritmo de Monte Carlo Hamiltoniano (HMC, Hamiltonian Monte Carlo) es una alternativa de alta eficiencia al algoritmo de Metropolis–Hastings.

El HMC incorpora una variable auxiliar denominada impulso (\(\boldsymbol{\varphi}\)), que permite recorrer la distribución objetivo mediante trayectorias dinámicas, superando las limitaciones del movimiento local propio de una caminata aleatoria.

La variable \(\boldsymbol{\varphi}\) incrementa la eficiencia del algoritmo al permitir desplazamientos más amplios y efectivos en el espacio de parámetros.

En HMC, las simulaciones se realizan sobre la distribución conjunta dada por \[ p(\boldsymbol{\theta}, \boldsymbol{\varphi} \mid \boldsymbol{y}) = p(\boldsymbol{\theta} \mid \boldsymbol{y}) \, p(\boldsymbol{\varphi} \mid \boldsymbol{y}) = p(\boldsymbol{\theta} \mid \boldsymbol{y}) \, p(\boldsymbol{\varphi}), \]
donde \(p(\boldsymbol{\varphi})\) suele ser una distribución Normal centrada. Aunque se simulan tanto \(\boldsymbol{\theta}\) como \(\boldsymbol{\varphi}\), el interés principal se centra en \(\boldsymbol{\theta}\).

Algoritmo

Inicialización del impulso:
Generar \(\boldsymbol{\varphi} \sim \mathsf{N}(\boldsymbol{0}, \mathbf{M})\), donde \(\mathbf{M}\) es la matriz de covarianza asociada a la función de impulso \(p(\boldsymbol{\varphi})\).
Integración Hamiltoniana:
Aplicar \(L\) pasos de integración con tamaño de paso \(\epsilon\), siguiendo el esquema leapfrog:
1. Primer medio paso en \(\boldsymbol{\varphi}\):
  \[ \boldsymbol{\varphi} \leftarrow \boldsymbol{\varphi} + \frac{\epsilon}{2} \, \frac{\partial}{\partial \boldsymbol{\theta}} \log p(\boldsymbol{\theta} \mid \boldsymbol{y}). \]
2. Paso completo en \(\boldsymbol{\theta}\)
  \[ \boldsymbol{\theta} \leftarrow \boldsymbol{\theta} + \epsilon \, \mathbf{M} \, \boldsymbol{\varphi}. \]
3. Pasos intermedios:
  Repetir la secuencia de actualización para \(\boldsymbol{\varphi}\) y \(\boldsymbol{\theta}\) durante \(L-1\) iteraciones. En cada paso intermedio:
  \[ \boldsymbol{\varphi} \leftarrow \boldsymbol{\varphi} + \epsilon \, \frac{\partial}{\partial \boldsymbol{\theta}} \log p(\boldsymbol{\theta} \mid \boldsymbol{y}), \]
  seguido de la actualización en \(\boldsymbol{\theta}\) como en (b).
4. Último medio paso en \(\boldsymbol{\varphi}\):
  \[ \boldsymbol{\varphi} \leftarrow \boldsymbol{\varphi} + \frac{\epsilon}{2} \, \frac{\partial}{\partial \boldsymbol{\theta}} \log p(\boldsymbol{\theta} \mid \boldsymbol{y}). \]
Cálculo de la probabilidad de aceptación:
Sean \((\boldsymbol{\theta}^{(b-1)}, \boldsymbol{\varphi}^{(b-1)})\) los valores iniciales y \((\boldsymbol{\theta}^*, \boldsymbol{\varphi}^*)\) los valores obtenidos tras \(L\) pasos. La razón de aceptación es:
\[ r = \frac{p(\boldsymbol{\theta}^* \mid \boldsymbol{y}) \, p(\boldsymbol{\varphi}^*)}{p(\boldsymbol{\theta}^{(b-1)} \mid \boldsymbol{y}) \, p(\boldsymbol{\varphi}^{(b-1)})}. \]
Aceptación o rechazo:
Actualizar:
\[ \boldsymbol{\theta}^{(b)} = \begin{cases} \boldsymbol{\theta}^*, & \text{con probabilidad } \min(1, r), \\ \boldsymbol{\theta}^{(b-1)}, & \text{en caso contrario}. \end{cases} \] Referencias:

MCMC using Hamiltonian dynamics: Disponible en arXiv. Este documento presenta los fundamentos teóricos del método y su aplicación en contextos de inferencia bayesiana.
Learning Hamiltonian Monte Carlo in R: Disponible en arXiv. Proporciona una introducción práctica al uso de HMC en R, con un enfoque en ejemplos y casos reales.
STAN Reference Manual: Disponible en STAN Documentation. Este manual describe en detalle la implementación de HMC en STAN, incluyendo configuraciones recomendadas y optimización de parámetros.

Consideraciones

Se recomienda seleccionar los valores de \(\mathbf{M}\), \(\epsilon\) y \(L\) de forma que la tasa de aceptación se sitúe entre el 60 % y el 70 %.
Una configuración habitual es \(\mathbf{M} = \mathbf{I}\), \(\epsilon = 1 / L\) y \(L\) como un entero positivo (e.g., \(L = 10\)).
No es necesario almacenar los valores de \(\boldsymbol{\varphi}\), ya que su función es únicamente facilitar un desplazamiento eficiente en el espacio de parámetros.

Ejemplo: Regresión Poisson (continuación)

Metropolis

# Hiperparámetros
beta0  <- rep(0, p)   # beta0 = 0
Sigma0 <- rep(10, p)  # Sigma0 = 10*I

# Parámetro de ajuste
Delta <- var(log(y + 1)) * solve(t(X)%*%X)

# Reproducibilidad
set.seed(123)

# Algoritmo de Metropolis
beta <- rep(0, p)                       # valor inicial
acr  <- 0                               # tasa de aceptación
B    <- 2100                            # número de muestras
BETA <- matrix(NA, nrow = B, ncol = p)  # almacenamiento

# Progreso
ncat <- floor(B / 10)

# Cadena
tictoc::tic()
for (b in 1:B) {
  # 1. Propuesta
  beta_p <- c(mvtnorm::rmvnorm(n = 1, mean = beta, sigma = Delta))

  # 2. Tasa de aceptación (escala log)
  eta_p <- as.numeric(X %*% beta_p)    # predictor lineal propuesto
  eta   <- as.numeric(X %*% beta)      # predictor lineal actual

  log_post_star <- sum(dpois(x = y, lambda = exp(eta_p), log = TRUE)) + 
       sum(dnorm(x = beta_p, mean = beta0, sd = Sigma0, log = TRUE))
       
  log_post_curr <- sum(dpois(x = y, lambda = exp(eta), log = TRUE)) + 
       sum(dnorm(x = beta, mean = beta0, sd = Sigma0, log = TRUE))

  log_r <- log_post_star - log_post_curr

  # 3) Actualizar
  if (runif(1) < exp(log_r)) {
    beta <- beta_p
    acr  <- acr + 1/B
  }

  # 4. Almacenar
  BETA[b, ] <- beta

  # 5. Progreso
  if (b %% ncat == 0) {
    cat(round(100 * b / B, 1), "% completado ...\n", sep = "")
  }
}

## 10% completado ...
## 20% completado ...
## 30% completado ...
## 40% completado ...
## 50% completado ...
## 60% completado ...
## 70% completado ...
## 80% completado ...
## 90% completado ...
## 100% completado ...

tictoc::toc()

## 0.151 sec elapsed

# Tasa de aceptación
round(100 * acr, 1)

## [1] 52.6

Monte Carlo Hamiltoniano

# Gradiente
grad <- function(beta, y, X) {
  out <- 0
  for (i in 1:length(y)) {
    eta <- sum(beta * X[i, ])
    out <- out + (y[i] - exp(eta)) * X[i, ]
  }
  return(out)
}

# Hiperparámetros
beta0  <- rep(0, p)    # beta0 = 0
Sigma0 <- diag(10, p)  # Sigma0 = 10*I

# Parámetros de ajuste
L   <- 100
eps <- 1/L

# Reproducibilidad
set.seed(123)

# Algoritmo de Metropolis
beta  <- rep(0, p)                       # valor inicial
acr   <- 0                               # tasa de aceptación
B     <- 2100                            # número de muestras
BETA2 <- matrix(NA, nrow = B, ncol = p)  # almacenamiento

# Cadena
tictoc::tic()
for (b in 1:B) {
  # 1. Inicialización del impulso
  phi    <- c(mvtnorm::rmvnorm(1, mean = rep(0, p), sigma = diag(p)))
  beta_p <- beta
  phi_p  <- phi
  
  # 2. Integración Hamiltoniana
  # Primer medio paso en phi
  phi_p <- phi_p + 0.5 * eps * grad(beta_p, y, X)
  
  # L pasos leapfrog
  for (l in 1:L) {
    # Paso completo en beta
    beta_p <- beta_p + eps * phi_p
    # Paso completo en phi (excepto al final)
    if (l < L) {
      phi_p <- phi_p + eps * grad(beta_p, y, X)
    }
  }
  
  # Último medio paso en phi
  phi_p <- phi_p + 0.5 * eps * grad(beta_p, y, X)
  
  # 3. Cálculo de la tasa de aceptación (log-scale)
  logr <- sum(dpois(y, lambda = exp(X %*% beta_p), log = TRUE)) -
          sum(dpois(y, lambda = exp(X %*% beta),   log = TRUE)) +
          mvtnorm::dmvnorm(beta_p, mean = beta0, sigma = Sigma0, log = TRUE) -
          mvtnorm::dmvnorm(beta,   mean = beta0, sigma = Sigma0, log = TRUE) +
          mvtnorm::dmvnorm(phi_p, log = TRUE) -
          mvtnorm::dmvnorm(phi,   log = TRUE)
  
  # 3. Aceptación o rechazo
  if (runif(1) < exp(logr)) {
    beta <- beta_p
    acr  <- acr + 1/B
  }
  
  # 4. Almacenamiento
  BETA2[b, ] <- beta
  
  # 5. Progreso
  if (b %% floor(B / 10) == 0) {
    cat(100 * round(b / B, 1), "% completado ...\n", sep = "")
  }
}

## 10% completado ...
## 20% completado ...
## 30% completado ...
## 40% completado ...
## 50% completado ...
## 60% completado ...
## 70% completado ...
## 80% completado ...
## 90% completado ...
## 100% completado ...

tictoc::toc()

## 9.704 sec elapsed

# Tasa de aceptación
round(100 * acr, 1)

## [1] 65

# Remover burn-in
burn <- 100
BETA  <- BETA[-seq_len(burn), , drop = FALSE]
BETA2 <- BETA2[-seq_len(burn), , drop = FALSE]

# Parámetros de gráfico
par(mfrow = c(3, 2), mar = c(2.75, 2.75, 1.5, 0.5), mgp = c(1.7, 0.7, 0))

for (j in seq_len(p)) {
  rango <- range(BETA[, j], BETA2[, j], na.rm = TRUE)

  # Metropolis
  plot(x = seq_len(nrow(BETA)), y = BETA[, j],
    col = 2, type = "p", pch = 16, cex = 0.6,
    ylim = rango, cex.axis = 0.9,
    main = "Metropolis", xlab = "Iteración",
    ylab = bquote(beta[.(j)]))

  # HMC
  plot(x = seq_len(nrow(BETA2)), y = BETA2[, j],
    col = 4, type = "p", pch = 16, cex = 0.6,
    ylim = rango, cex.axis = 0.9,
    main = "Hamiltoniano", xlab = "Iteración",
    ylab = bquote(beta[.(j)]))
}

# Tamaños efectivos (neff) y tabla comparativa
neff_metro <- round(coda::effectiveSize(BETA),  1)  # vector por parámetro
neff_hmc   <- round(coda::effectiveSize(BETA2), 1)

# Armar tabla
tab <- cbind("neff Metro" = neff_metro, "neff HMC" = neff_hmc)

# Nombres de las filas
rownames(tab) <- paste0("beta", seq_len(p))

# Mostrar
round(tab, 3)

##       neff Metro neff HMC
## beta1      118.5   1176.3
## beta2      110.9   1154.9
## beta3      103.9   1137.9

Ejercicios

Regresión logística: Ejercicio 10.2 de Hoff (2009).
Regresión con errores correlacionados: Ejercicio 10.2 de Hoff (2009).
Selección de modelos: Ejercicio 10.5 de Hoff (2009).

Referencias

Hoff, P. D. (2009). A First Course in Bayesian Statistical Methods. Springer New York.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). Chapman & Hall/CRC.

Regresión Poisson

Juan Sosa PhD

Email jcsosam@unal.edu.co

GitHub https://github.com/jstats1702

Modelo

Distibución conjunta posterior

Algorithmo de Metropolis-Hastings

Algoritmo de Metropolis

Algoritmo

Consideraciones

Ejemplo: Modelo Normal

Algoritmo de Metropolis adaptativo

Algoritmo Metropolis-Hastings

Algoritmo

Consideraciones

Ejemplo: Regresión Poisson

Tratamiento de datos

Análisis exploratorio

Ajuste del modelo

Algoritmo de Monte Carlo Hamiltoniano

Algoritmo

Consideraciones

Ejemplo: Regresión Poisson (continuación)

Metropolis

Monte Carlo Hamiltoniano

Ejercicios

Referencias