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).
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Ć:
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.
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.
\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\,.
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.
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.
Dado el estado actual \theta^{(b)} del parƔmetro \theta, el siguiente estado \theta^{(b+1)} se genera como sigue:
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.
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).
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}
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.
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 de propuestas
arbitraria.
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.
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}
Algoritmo de Metropolis en un modelo Normal con varianza conocida.
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.
Se considera una distribuciĆ³n previa Normal para el parĆ”metro \theta: \theta \sim \textsf{Normal}(\mu, \tau^2), donde:
# 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
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
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.
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}).
ActualizaciĆ³n mediante integraciĆ³n Hamiltoniana: Aplicar L pasos de integraciĆ³n escalados por un tamaƱo de paso \epsilon:
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}).
Paso completo en \boldsymbol{\theta}: \boldsymbol{\theta} \leftarrow \boldsymbol{\theta} + \epsilon \, \mathbf{M} \, \boldsymbol{\varphi}.
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}.
Ć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 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)})}.
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\}.
# 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:
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.
# 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