Processing math: 1%

1 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

yiāˆ£Īøiiidāˆ¼Poisson(Īøi),i=1,ā€¦,n, donde la media Īø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}^k \beta_j x_{i,j}.

AquĆ­:

  • \boldsymbol{\beta} = (\beta_1, \ldots, \beta_k): vector de coeficientes del modelo.
  • \boldsymbol{x}_i = (x_{i,1}, \ldots, x_{i,k}): vector de predictores para la i-Ć©sima observaciĆ³n.

En este modelo, la funciĆ³n de enlace (link function) es la funciĆ³n logarĆ­tmica (\log), que garantiza que \theta_i > 0, cumpliendo con las propiedades de la distribuciĆ³n Poisson.

DistribuciĆ³n Previa

Dado que las distribuciones previas conjugadas para modelos lineales generalizados no existen (excepto para el modelo de regresiĆ³n Normal), se utiliza una distribuciĆ³n previa Normal como aproximaciĆ³n:

\boldsymbol{\beta} \sim \textsf{N}(\boldsymbol{\beta}_0, \Sigma_0), donde: - \boldsymbol{\beta}_0: Media previa de los coeficientes. - \Sigma_0: Matriz de covarianza previa.

Componentes del Modelo

  1. ParƔmetros del modelo: (\beta_1, \ldots, \beta_k), que representan los efectos de los predictores sobre \log(\theta_i).
  2. HiperparĆ”metros del modelo: (\boldsymbol{\beta}_0, \Sigma_0), que especifican las propiedades de la distribuciĆ³n previa de \boldsymbol{\beta}.

2 DistibuciĆ³n conjunta posterior

\begin{align*} p(\boldsymbol{\beta}\mid\boldsymbol{y}) &\propto p(\boldsymbol{y}\mid\boldsymbol{\beta})\,p(\boldsymbol{\beta}) \\ &= \prod_{i=1}^n \textsf{Poisson}(y_i\mid\theta_i) \times\textsf{N}(\boldsymbol{\beta}\mid\boldsymbol{\beta}_0,\Sigma_0) \\ &= \prod_{i=1}^n \frac{e^{-\theta_i}\,\theta_i^{y_i}}{y_i!} \times (2\pi)^{-k/2} \, |\Sigma_0|^{-1/2} \, \exp\left\{ -\tfrac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_0)^{\textsf{T}}\Sigma_0^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_0) \right\} \\ &\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\} \end{align*} con \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.

Por lo tanto, 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}).

AdemƔs, 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\,.

3 Algorithmo de Metropolis-Hastings

El algoritmo de Metropolis-Hastings es un mĆ©todo general para aproximar distribuciones posteriores de cualquier forma, incluso cuando no se dispone de una soluciĆ³n analĆ­tica exacta.

El desafĆ­o principal surge cuando p(\boldsymbol{\theta} \mid \boldsymbol{y}) no corresponde a una distribuciĆ³n estĆ”ndar conocida, lo que imposibilita generar muestras directamente de ella de manera sencilla.

3.1 Algoritmo de Metropolis

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

Este algoritmo genera una cadena de Markov que converge hacia la distribuciĆ³n posterior deseada p(\theta \mid \boldsymbol{y}), asumiendo 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 como sigue:

  1. 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 la distribuciĆ³n de propuestas (proposal distribution), suele ser: 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, seleccionado para garantizar un buen rendimiento del algoritmo.

  2. CĆ”lculo de la tasa de aceptaciĆ³n:
    Calcular la tasa de aceptaciĆ³n r, definida como: r = \frac{p(\theta^* \mid \boldsymbol{y})}{p(\theta^{(b)} \mid \boldsymbol{y})}. Para evitar problemas de 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).

  3. ActualizaciĆ³n del estado:
    Establecer el siguiente 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

  • La simetrĆ­a de J(\theta \mid \theta^{(b)}) simplifica el cĆ”lculo de la tasa de aceptaciĆ³n, ya que no es necesario considerar las probabilidades de la distribuciĆ³n de propuestas en el cĆ”lculo, lo que reduce la complejidad computacional.
  • 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.
  • Este algoritmo genera 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.
  • 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.
  • Es comĆŗn calibrar \delta para lograr una tasa de aceptaciĆ³n entre el 30% y el 50%, lo que suele optimizar el rendimiento del algoritmo. Para facilitar esta calibraciĆ³n, se pueden utilizar algoritmos adaptativos que ajustan automĆ”ticamente el valor de \delta durante las iteraciones iniciales.
  • El algoritmo de Metropolis es particularmente Ćŗtil para simular distribuciones posteriores complejas que no tienen formas analĆ­ticas conocidas, proporcionando una herramienta robusta para la inferencia en modelos Bayesianos.

3.2 Algoritmo Metropolis-Hastings

El algoritmo Metropolis-Hastings es una generalizaciĆ³n de los algoritmos de Gibbs y 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 genera de la siguiente manera:

  1. 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.

  2. CĆ”lculo de la tasa de aceptaciĆ³n:
    Calcular 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 cĆ”lculo ajusta la probabilidad de aceptaciĆ³n para tener en cuenta la no simetrĆ­a de la distribuciĆ³n de propuestas J.

  3. ActualizaciĆ³n del estado:
    Actualizar el estado \theta^{(b+1)} segĆŗ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 especial de Metropolis-Hastings cuando 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).
  • El algoritmo Metropolis-Hastings permite abordar problemas en los que las distribuciones condicionales completas del Gibbs ni las propuestas simĆ©tricas del Metropolis son adecuadas.

4 Ejemplo: Modelo Normal

Algoritmo de Metropolis en un modelo Normal con varianza conocida.

DistribuciĆ³n muestral

Los datos y_i se asumen independientes e idĆ©nticamente distribuidos (i.i.d.) segĆŗn una distribuciĆ³n Normal con media \theta y varianza 1: y_i \mid \theta \stackrel{\text{iid}}{\sim} \textsf{Normal}(\theta, 1), \qquad i = 1, \ldots, n.

DistribuciĆ³n previa

Se considera una distribuciĆ³n previa Normal para el parĆ”metro \theta: \theta \sim \textsf{Normal}(\mu, \tau^2), donde:

  • \mu es la media previa.
  • \tau^2 es la varianza previa, que refleja la incertidumbre inicial sobre \theta.
# simular datos
n  <- 5
s2 <- 1 
set.seed(123)
y <- round(rnorm(n = n, mean = 10, sd = 1), 2)
# previa
t2 <- 10 
mu <- 5
# Metropolis
theta <- 0       # valor inicial
delta <- 1.75    # parƔmetro de ajuste
S     <- 100000  # numero de muestras
THETA <- NULL    # almacenamiento
# cadena
set.seed(123)
for (s in 1:S) {
  # 1. propuesta
  theta.star <- rnorm(1,theta,sqrt(delta))
  # 2. tasa de aceptaciĆ³n
  log.r <- (sum(dnorm(y, theta.star, sqrt(s2), log = T)) + dnorm(theta.star, mu, sqrt(t2), log = T)) - (sum(dnorm(y, theta, sqrt(s2), log = T)) + dnorm(theta, mu, sqrt(t2), log = T)) 
  # 3. actualizar
  if (runif(1) < exp(log.r)) 
    theta <- theta.star 
  # 4. almacenar
  THETA <- c(THETA, theta)
}

ACR    <- NULL  # tasa de aceptaciones
ACF    <- NULL  # autocorrelaciones
THETAA <- NULL  # muestras
for(delta2 in 2^c(-5,-1,1,5,7) ) {
  # parƔmetros iniciales
  THETA <- NULL
  S     <- 100000
  theta <- 0
  acs   <- 0  # tasa de aceptaciĆ³n
  # cadena
  set.seed(1)
  for(s in 1:S) {
    # 1. propuesta
    theta.star <- rnorm(1,theta,sqrt(delta2))
    # 2. tasa de aceptaciĆ³n
    log.r <- (sum(dnorm(y, theta.star, sqrt(s2), log = T)) + dnorm(theta.star, mu, sqrt(t2), log = T)) - (sum(dnorm(y, theta, sqrt(s2), log = T)) + dnorm(theta, mu, sqrt(t2), log = T)) 
    # 3. actualizar
    if(runif(1) < exp(log.r)) { 
      theta <- theta.star 
      acs   <- acs + 1 
    }
    # 4. almacenar
    THETA <- c(THETA, theta) 
  }
  # almacenar valores de todos los casos (delta2)
  ACR    <- c(ACR, acs/s) 
  ACF    <- c(ACF, acf(THETA, plot = F)$acf[2])
  THETAA <- cbind(THETAA, THETA)
}

# tasas de aceptaciĆ³n
round(ACR, 3)
## [1] 0.874 0.572 0.357 0.098 0.050
# autocorrelaciones
round(ACF ,3)
## [1] 0.950 0.683 0.649 0.870 0.929
# tamaƱos efectivos de muestra
round(coda::effectiveSize(THETAA), 3)
##     THETA     THETA     THETA     THETA     THETA 
##  2486.644 15461.739 19934.233  6466.263  3334.037

5 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.

Referencia: 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.

Las variables del modelo son:

# Actividades de reproducciĆ³n de gorriones (hembras) en funciĆ³n de la edad
# "age"     : edad.
# "fledged" : nĆŗmero de crĆ­as.
spfage <- 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"))
)
spfage <- as.data.frame(spfage)
# diseƱo
spf  <- spfage$fledged  # y  = variable dependiente (respuesta)
age  <- spfage$age      # x1 = variable independiente 1
age2 <- age^2           # x2 = variable independiente 2
# visualizaciĆ³n
par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot(spf~as.factor(age), range=0, xlab="Edad (aƱos)", ylab="No. Crias", col="gray", border="lightgray")

# GLM (frecuentista)
summary(glm(spf~age+age2,family="poisson"))
## 
## Call:
## glm(formula = spf ~ 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

En muchos problemas, utilizar la varianza posterior como distribuciĆ³n de propuestas es una elecciĆ³n eficiente. Aunque la varianza posterior exacta no se conoce antes de ejecutar el algoritmo de Metropolis, una aproximaciĆ³n razonable suele ser suficiente.

Si la tasa de aceptaciĆ³n resultante es demasiado alta o demasiado baja, se puede ajustar la variabilidad de la propuesta de manera iterativa para mejorar el rendimiento del algoritmo. Este enfoque permite equilibrar la exploraciĆ³n eficiente del espacio de parĆ”metros con la convergencia adecuada.

# datos 
y <- spf                                # variable respuesta
X <- cbind(rep(1,length(y)),age,age^2)  # matriz de diseƱo
n <- dim(X)[1]                          # tamaƱo de la muestra
p <- dim(X)[2]                          # numero de predictores
# previa
pmn.beta <- rep(0, p)         # beta0 = 0
psd.beta <- rep(sqrt(10), p)  # Sigma0 = 10*I
# matriz de varianza para la propuesta
# y + 1/2 evitar problemas en la frontera con 0
var.prop <- var(log(y + 1)) * solve(t(X)%*%X)
# algoritmo de Metropolis
S    <- 100000
BETA <- matrix(data = NA, nrow = S, ncol = p)
ac   <- 0
ncat <- floor(S/10)
beta <- rep(0, p) 
# cadena
set.seed(123)
for(s in 1:S) {
  # 1. propuesta
  beta.p <- c(mvtnorm::rmvnorm(n = 1, mean = beta, sigma = var.prop))
  # 2. tasa de aceptaciĆ³n
  lhr <- sum(dpois(x = y, lambda = exp(X%*%beta.p), log = T)) - sum(dpois(x = y, lambda = exp(X%*%beta), log = T)) + sum(dnorm(x = beta.p, mean = pmn.beta, sd = psd.beta, log = T)) - sum(dnorm(x = beta, mean = pmn.beta, sd = psd.beta, log = T))
  # 3. actualizar
  if (runif(1) < exp(lhr)) { 
    beta <- beta.p 
    ac   <- ac + 1 
  }
  # 4. almacenar
  BETA[s,] <- beta
  # 5. Progreso
  if (s%%ncat == 0)
    cat(100*round(s/S, 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*ac/S, 1)
## [1] 52.6
# tamaƱo efectivo de muestra
round(apply(X = BETA, MARGIN = 2, FUN = coda::effectiveSize), 1)
## [1] 7044.6 6065.5 5522.3

round(quantile(BETA[,2], c(.025, .975)), 3)
##  2.5% 97.5% 
## 0.080 1.388
round(quantile(BETA[,3], c(.025, .975)), 3)
##   2.5%  97.5% 
## -0.257 -0.032
round(mean(BETA[,2] > 0), 3)
## [1] 0.986
round(mean(BETA[,3] > 0), 3)
## [1] 0.005

6 Algoritmo de Monte Carlo Hamiltoniano

El algoritmo de Monte Carlo Hamiltoniano (HMC, Hamiltonian Monte Carlo) es una alternativa altamente eficiente al Algoritmo de Metropolis-Hastings. Este mĆ©todo introduce una variable auxiliar, denominada ā€œimpulsoā€ (\boldsymbol{\varphi}), que permite explorar la distribuciĆ³n objetivo de manera mĆ”s efectiva mediante trayectorias dinĆ”micas, evitando las limitaciones del movimiento de caminata aleatoria local inherente a otros mĆ©todos.

En HMC, las simulaciones se realizan sobre la distribuciĆ³n conjunta: p(\boldsymbol{\theta}, \boldsymbol{\varphi} \mid \boldsymbol{y}) = p(\boldsymbol{\theta} \mid \boldsymbol{y}) \, p(\boldsymbol{\varphi}), donde el tĆ©rmino p(\boldsymbol{\varphi}) corresponde a una distribuciĆ³n normal centrada. Aunque ambas variables, \boldsymbol{\theta} y \boldsymbol{\varphi}, son simuladas, el interĆ©s principal recae Ćŗnicamente en \boldsymbol{\theta}. La variable \boldsymbol{\varphi} actĆŗa como un auxiliar que mejora la eficiencia del algoritmo al facilitar desplazamientos mĆ”s amplios en el espacio de parĆ”metros.

Algoritmo HMC

  1. InicializaciĆ³n del impulso: Simular \boldsymbol{\varphi} \sim \text{N}(\boldsymbol{0}, \mathbf{M}), donde \mathbf{M} es una matriz diagonal que define la covarianza de la funciĆ³n de impulso p(\boldsymbol{\varphi}).

  2. ActualizaciĆ³n mediante integraciĆ³n Hamiltoniana: Aplicar L pasos de integraciĆ³n escalados por un tamaƱo de paso \epsilon:

    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. Siguientes pasos alternados: Repetir el proceso de actualizaciĆ³n de \boldsymbol{\varphi} y \boldsymbol{\theta} L-1 veces. En cada iteraciĆ³n: \boldsymbol{\varphi} \leftarrow \boldsymbol{\varphi} + \epsilon \, \frac{\partial}{\partial \boldsymbol{\theta}} \log p(\boldsymbol{\theta} \mid \boldsymbol{y}). En la Ćŗltima iteraciĆ³n, solo es necesario completar el paso en \boldsymbol{\theta}.

    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}).

  3. CĆ”lculo de la tasa de aceptaciĆ³n: Sean (\boldsymbol{\theta}^{(b-1)}, \boldsymbol{\varphi}^{(b-1)}) los valores iniciales y (\boldsymbol{\theta}^*, \boldsymbol{\varphi}^*) los valores despuĆ©s de los L pasos. Calcular: r = \frac{p(\boldsymbol{\theta}^* \mid \boldsymbol{y}) \, p(\boldsymbol{\varphi}^*)}{p(\boldsymbol{\theta}^{(b-1)} \mid \boldsymbol{y}) \, p(\boldsymbol{\varphi}^{(b-1)})}.

  4. AceptaciĆ³n o rechazo: Establecer \boldsymbol{\theta}^{(b)} = \boldsymbol{\theta}^* con probabilidad \min\{r, 1\}, o \boldsymbol{\theta}^{(b)} = \boldsymbol{\theta}^{(b-1)} con probabilidad 1 - \min\{r, 1\}.

Consideraciones

  • Se recomienda elegir los valores de \mathbf{M}, \epsilon y L para obtener una tasa de aceptaciĆ³n entre el 60% y 70%.
  • Una configuraciĆ³n comĆŗn es \mathbf{M} = \mathbf{I}, \epsilon = 1 / L y L como un entero positivo (por ejemplo, L = 10).
  • No es necesario almacenar los valores de \boldsymbol{\varphi}, ya que su Ćŗnico propĆ³sito es facilitar el desplazamiento eficiente en el espacio de parĆ”metros.

7 Ejemplo: RegresiĆ³n Poisson (continuaciĆ³n)

# datos 
y <- spf                                # variable respuesta
X <- cbind(rep(1,length(y)),age,age^2)  # matriz de diseƱo
n <- dim(X)[1]                          # tamaƱo de la muestra
p <- dim(X)[2]                          # numero de predictores
# previa
pmn.beta <- rep(0, p)         # beta0 = 0
psd.beta <- rep(sqrt(10), p)  # Sigma0 = 10*I
# matriz de varianza para la propuesta
# y + 1/2 evitar problemas en la frontera con 0
var.prop <- var(log(y + 1)) * solve(t(X)%*%X)
# algoritmo de Metropolis
S    <- 2100
BETA <- matrix(data = NA, nrow = S, ncol = p)
ac   <- 0
ncat <- floor(S/10)
beta <- rep(0, p) 
# cadena
set.seed(123)
for(s in 1:S) {
  # 1. propuesta
  beta.p <- c(mvtnorm::rmvnorm(n = 1, mean = beta, sigma = var.prop))
  # 2. tasa de aceptaciĆ³n
  lhr <- sum(dpois(x = y, lambda = exp(X%*%beta.p), log = T)) - sum(dpois(x = y, lambda = exp(X%*%beta), log = T)) + sum(dnorm(x = beta.p, mean = pmn.beta, sd = psd.beta, log = T)) - sum(dnorm(x = beta, mean = pmn.beta, sd = psd.beta, log = T))
  # 3. actualizar
  if (runif(1) < exp(lhr)) { 
    beta <- beta.p 
    ac   <- ac + 1 
  }
  # 4. almacenar
  BETA[s,] <- beta
  # 5. Progreso
  if (s%%ncat == 0)
    cat(100*round(s/S, 1), "% completado ... \n", sep = "" )
}
## 10% completado ... 
## 20% completado ... 
## 30% completado ... 
## 40% completado ... 
## 50% completado ... 
## 60% completado ... 
## 70% completado ... 
## 80% completado ... 
## 90% completado ... 
## 100% completado ...

Hamiltonian Monte Carlo (HMC)

Para una explicaciĆ³n detallada y ejemplos prĆ”cticos de este mĆ©todo, se pueden consultar las siguientes referencias:

# gradiente
grad <- function(bet, y, X) {
  out <- 0
  for (i in 1:length(y))
    out <- out + (y[i] - exp(sum(bet*X[i,])))*X[i,]
  out
}
# previa
mu0    <- rep(0,  p)   # beta0 = 0
Sigma0 <- diag(10, p)  # Sigma0 = 10*I
# cadena
L   <- 100
eps <- 1/L
bet <- rep(0, p) # valor inicial beta
ac  <- 0         # tasa de aceptaciĆ³n
S   <- 2100
BETA2 <- matrix(data = NA, nrow = S, ncol = p)
set.seed(123)
for(s in 1:S) {
  # 1. propuesta
  bet.p <- bet
  z <- c(mvtnorm::rmvnorm(n = 1, mean = rep(0,p), sigma = diag(p)))
  z.p <- z
  z.p <- z.p + 0.5*eps*grad(bet.p, y, X)
  for (l in 1:L) {
    bet.p <- bet.p + eps*z.p
    if (l < L)
      z.p   <- z.p + eps*grad(bet.p, y, X)
  }
  z.p <- z.p + 0.5*eps*grad(bet.p, y, X)
  # 2. tasa de aceptaciĆ³n
  logr <- sum(dpois(x = y, lambda = exp(X%*%bet.p), log = T)) - sum(dpois(x = y, lambda = exp(X%*%bet),log = T)) + mvtnorm::dmvnorm(x = bet.p, mean = mu0, sigma = Sigma0, log = T) - mvtnorm::dmvnorm(x = bet, mean = mu0, sigma = Sigma0, log = T) + mvtnorm::dmvnorm(x = z.p, log = T) - mvtnorm::dmvnorm(x = z, log = T)
  # 3. actualizar
  if (runif(n = 1) < exp(logr)) { 
    bet <- bet.p 
    ac  <- ac + 1 
  }
  # 4. almacenar
  BETA2[s,] <- bet
  # 5. Progreso
  if (s%%floor(S/10) == 0) cat(100*round(s/S, 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(ac/S, 2)
## [1] 0.65
BETA <- BETA[-c(1:100),]
BETA2 <- BETA2[-c(1:100),]
par(mfrow = c(3,2), mar = c(2.75,2.75,1.5,.5), mgp = c(1.7,.7,0))
for (j in 1:p) {
  rango <- range(BETA[,j], BETA2[,j])
  # Metropolis
  plot(x = 1:nrow(BETA), y = BETA[,j], col = 2, type = "p", pch = 16, ylim = rango, cex.axis = 0.9, main = "Metropolis", xlab = "IteraciĆ³n", ylab = "")
  # HMC
  plot(x = 1:nrow(BETA2), y = BETA2[,j], col = 4, type = "p", pch = 16, ylim = rango, cex.axis = 0.9, main = "Hamiltoniano", xlab = "IteraciĆ³n", ylab = "")
}

# tamaƱos efectivos & CV de MC
tab <- cbind(
  round(apply(X = BETA,  MARGIN = 2, FUN = coda::effectiveSize), 1),
  round(apply(X = BETA2, MARGIN = 2, FUN = coda::effectiveSize), 1)
)
rownames(tab) <- paste("beta", 1:3)
colnames(tab) <- c("neff Metro", "neff HMC")
round(tab, 3)
##        neff Metro neff HMC
## beta 1      123.3   1176.3
## beta 2      117.9   1154.9
## beta 3      110.5   1137.9

Referencias