MCMC
jmm
23 de abril de 2019
Vamos a hacer un ejemplo simple de MCMC para ajustar una regresión.
Primerlo simulamos datos
set.seed(1)
n <- 30
b0 <- 0
b1 <- 0.9
x <- runif(n, -2, 2)
y <- rpois(n, lambda = exp(b0 + b1 * x))
plot(x, y)
Vamos a tener que calcular muchas veces la probabilidad de los datos condicional en los parámetros (\(p(y| \theta)\)), pero para evitar problemas numéricos trabajamos con el logaritmo de esa probabilidad, así que definimos una función para calcular el logaritmo de likelihood de una combinación de parámetros. Aquí es importante acordarse de usar log = TRUE.
like <- function(x, y, b0, b1){
lambda <- exp(b0 + b1 * x)
ll <- sum(dpois(y, lambda = lambda, log = TRUE))
}Definimos cuántas iteraciones de MCMC queremos usar, Asumimos previas con distribución normal para los coeficientes de la regresión y definimos los valores de media y desvío para estas previas.
Vamos a usar un MCMC tipo “random walk”, entonces tenemos que definir el desvió de la distribución de propuestas
n.iter <- 10000
m.b0 <- 0
m.b1 <- 0
sd.b0 <- 5
sd.b1 <- 1
sd.jump <- 1
sims <- matrix(NA, n.iter, 2)
pos <- numeric(n.iter)Ahora definimos valores iniciales
new.b0 <- runif(1, -1, 1)
new.b1 <- runif(1, -1, 1)
pos[1] <- like(x, y, new.b0, new.b1) +
dnorm(new.b0, m.b0, sd.b0, log = TRUE) +
dnorm(new.b1, m.b1, sd.b1, log = TRUE)
sims[1, 1] <- new.b0
sims[1, 2] <- new.b1y ahora el MCMC
for(i in 2:n.iter){
# si no aceptamos nuevos valores, todo queda como antes
pos[i] <- pos[i-1]
sims[i,] <- sims[i-1,]
# update b0
new.b0 <- rnorm(1, sims[i-1, 1], sd.jump)
n.pos <- like(x, y, new.b0, sims[i-1, 2]) +
dnorm(new.b0, m.b0, sd.b0, log = TRUE) +
dnorm(sims[i-1, 2], m.b1, sd.b1, log = TRUE)
if( (n.pos - pos[i]) > log(runif(1)) ){
pos[i] <- n.pos
sims[i,1] <- new.b0
}
#update b1
new.b1 <- rnorm(1, sims[i-1,2], sd.jump)
n.pos <- like(x, y, sims[i, 1], new.b1) +
dnorm(sims[i, 1], m.b0, sd.b0, log = TRUE) +
dnorm(new.b1, m.b1, sd.b1, log = TRUE)
if( (n.pos - pos[i]) > log(runif(1)) ){
pos[i] <- n.pos
sims[i, 2] <- new.b1
}
}Ahora graficamos las cadenas luego de descartar el “burn in”
burn = ceiling(n.iter/2)
op <- par(mfrow = c(1,2), las = 1, cex.lab = 1.3)
plot(sims[burn:n.iter, 1], type = "l", ylab = expression(b[0]))
abline(h = b0)
plot(sims[burn:n.iter, 2], type = "l", ylab = expression(b[1]))
abline(h = b1)
par(op)Ahora la posterior conjunta
library(ggplot2)
library(ggExtra)
#library(gridExtra)
df <- data.frame(b0 = sims[burn:n.iter,1], b1 = sims[burn:n.iter,2])
p <- ggplot(df, aes(b0, b1)) + geom_point(alpha = 0.1) + theme_classic() + xlim(-1, 1) + ylim(0, 1.5) + coord_equal()
ggExtra::ggMarginal(p, type = "histogram")## Warning: Removed 5 rows containing missing values (geom_point).
Ejercicios
1- Cómo cambian las posteriores cuando cambia el número de datos? Probar con n = 5, n = 100 y n = 500.
2- Modifique el scrip anterior para contar cuántas veces se acepta mover la cadena Markpviana para cada parámetro
3- Qué valor de sd.jump resulta en aproximadamente 30% de aceptación en movimientos de las cadenas?
Ejemplo con más de una cadena
n.iter <- 10000
n.chain <- 3
m.b0 <- 0
m.b1 <- 0
sd.b0 <- 1
sd.b1 <- 1
sd.jump <- 1
pos <- matrix(NA, n.iter, n.chain)
sims <- vector("list", n.chain)
for(i in 1: n.chain){
sims[[i]] <- matrix(NA, n.iter, 2)
}Definimos valores iniciales
for(i in 1:n.chain){
new.b0 <- runif(1, -1, 1)
new.b1 <- runif(1, -1, 1)
pos[1,i] <- like(x, y, new.b0, new.b1) +
dnorm(new.b0, m.b0, sd.b0) +
dnorm(new.b1, m.b1, sd.b1)
sims[[i]][1, 1] <- new.b0
sims[[i]][1, 2] <- new.b1
}Y corremos las cadenas
for(i in 2:n.iter){
for(j in 1: n.chain){
# si no aceptamos nuevos valores, todo queda como antes
pos[i,j] <- pos[i-1,j]
sims[[j]] [i,] <- sims[[j]] [i-1,]
# update b0
new.b0 <- rnorm(1, sims[[j]] [i-1, 1], sd.jump)
n.pos <- like(x, y, new.b0, sims[[j]] [i-1, 2]) +
dnorm(new.b0, m.b0, sd.b0) +
dnorm(sims[[j]] [i-1, 2], m.b1, sd.b1)
if( (n.pos - pos[i,j]) > log(runif(1)) ){
pos[i,j] <- n.pos
sims[[j]] [i,1] <- new.b0
}
#update b1
new.b1 <- rnorm(1, sims[[j]] [i-1,2], sd.jump)
n.pos <- like(x, y, sims[[j]] [i, 1], new.b1) +
dnorm(sims[[j]] [i, 1], m.b0, sd.b0) +
dnorm(new.b1, m.b1, sd.b1)
if( (n.pos - pos[i,j]) > log(runif(1)) ){
pos[i,j] <- n.pos
sims[[j]] [i, 2] <- new.b1
}
}
}Ahora usamos el paquete coda para ver si hay convergencia
library(coda)
burn <- ceiling(n.iter/2)
thin <- 1
sim <- vector("list", n.chain)
ps <- mcmc(pos, start = burn, end = n.iter, thin = thin)
for(i in 1:n.chain){
sim[[i]] <- as.mcmc(cbind(mcmc(sims[[i]], start = burn, end = n.iter, thin = thin),
ps[,i]))
colnames(sim[[i]]) <- list("b0", "b1", "pos")
}
gelman.diag(as.mcmc.list(sim))## Potential scale reduction factors:
##
## Point est. Upper C.I.
## b0 1.02 1.08
## b1 1.02 1.07
## pos 1.01 1.02
##
## Multivariate psrf
##
## 1.02También podemos ver el tamaño efectivo de las muestras de las posteriores
effectiveSize(sim)## b0 b1 pos
## 420.6666 388.8462 1143.6892La tasa de rechazo de las propuestas de Metropolis-Hastings
rejectionRate(as.mcmc.list(sim))## b0 b1 pos
## 0.8234667 0.8694000 0.7159333y usar las gráficas de coda
plot(as.mcmc(sim))
Para ver las medias de las posteriores y los cuantiles de las posteriores:
summary(as.mcmc(sim))##
## Iterations = 1:5001
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 5001
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## b0 -0.1046 0.2311 0.001887 0.010765
## b1 0.9339 0.1623 0.001325 0.008092
## pos -36.7761 1.1645 0.009507 0.034328
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## b0 -0.6052 -0.2519 -0.08497 0.06058 0.2926
## b1 0.6470 0.8185 0.92301 1.04097 1.2920
## pos -39.3665 -37.0959 -36.48287 -36.08987 -35.8121Para calcular los intervalos de credibilidad usando “highest posterior density”:
S <- as.mcmc(rbind(sim[[1]], sim[[2]], sim[[3]]))
HPDinterval(S)## lower upper
## b0 -0.5863218 0.2927372
## b1 0.6239169 1.2375000
## pos -38.7353574 -35.7782172
## attr(,"Probability")
## [1] 0.95001Just Another Gibbs Sampler
JAGS es un software que programa cadenas de Markov Chain Monte Carlo (MCMC) para modelos bayesianos Plummer. JAGS es un sucesor de BUGS, que es [Bayesian inference using Gibbs sampling] (https://www.mrc-bsu.cam.ac.uk/software/bugs/). JAGS es muy parecido a BUGS pero tiene algunas funciones extra y a veces es más rápido. Además, JAGS se puede usar en Windows, Mac y Linux.
Para usar JAGS (o BUGS) primero tenemos que definir el modelo usando un lenguaje sintético y guardarlo como un archivo de texto. Podemos escribir el modelo en un editor de texto o, como en el ejemplo de abajo, usar la función cat de R que nos permite escribir y guardar texto desde un script.
cat(file = "rp.bug",
"
model{
# likelihood
for (i in 1:n) {
y[i] ~ dpois(lambda[i])
log(lambda[i]) <- b0 + b1 * x[i]
}
# previas
b0 ~ dnorm(0, 0.04)
b1 ~ dnorm(0, 1)
}"
)Además, tenemos que armar una lista con los datos que vamos a pasarle a JAGS, una función para generar los valores iniciales de las cadenas Markovianas, y definir los parámetros que queremos guardar.
m.data <- list("n", "y", "x")
inits <- function() list(b0 = runif(1,-1,1),
b1 = runif(1, -1,1))
params <- c("b0", "b1")Finalmente, definimos cuántas iteraciones queremos ejecutar por cada cadena, cuántos valores vamos a ir descartado en una secuencia (“thin”), qué largo tiene el “burn in” y cuántas cadenas queremos correr.
ni <- 10000 # número de iteraciones
nt <- 1 # tasa de "thining"
nb <- 5000 # cuantas iteraciones usamos de "burn in"
nc <- 5 # y cuantas cadenas corremosPara llamar a JAGS desde R usamos el paquete jagsUI y la función jags de ese paquete
library(jagsUI)## Loading required package: lattice##
## Attaching package: 'jagsUI'## The following object is masked from 'package:coda':
##
## traceplot## The following object is masked from 'package:utils':
##
## Viewjags.sim <- jags(data = m.data,
inits = inits,
parameters.to.save = params,
model.file = "rp.bug",
n.chains = nc,
n.thin = nt,
n.iter = ni,
n.burnin = nb)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 30
## Unobserved stochastic nodes: 2
## Total graph size: 157
##
## Initializing model
##
## Adaptive phase.....
## Adaptive phase complete
##
##
## Burn-in phase, 5000 iterations x 5 chains
##
##
## Sampling from joint posterior, 5000 iterations x 5 chains
##
##
## Calculating statistics.......
##
## Done.Para ver los resultados de este análisis podemos pedirle a R que “imprima” la salida de la función jags usando print(m.sim). Vemos que se reporta el nombre del modelo que se usó, cuántas cadenas se simularon y otros detalles como el tiempo de ejecución. Después aparece una tabla con la lista de parámetros que le pedimos que registre y la devianza. En la tabla aparece la media, desvío y cuantiles estimados a partir de las cadenas Markovianas. También aparece una columna (overlap0) que nos dice si la posterior incluye al cero o no, y otra (f) que nos dice qué fracción de la posterior es del mismo signo que la media. Finalmente, aparecen dos columnas con información importante. Rhat estima si las cadenas convergieron a una distribución estable y n.eff estima el número efectivo de muestras de la posterior que surgen de las cadenas. Antes de sacar cualquier conclusión a partir de la salida de JAGS, tenemos que chequear que las cadenas hayan convergido (Rhat \(\leq 1.1\)). También, ver que el tamaño efectivo de la muestra de la posterior (n.eff) sea suficiente.
print(jags.sim)## JAGS output for model 'rp.bug', generated by jagsUI.
## Estimates based on 5 chains of 10000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 5000 iterations and thin rate = 1,
## yielding 25000 total samples from the joint posterior.
## MCMC ran for 0.049 minutes at time 2019-04-26 12:26:42.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## b0 -0.104 0.230 -0.584 -0.093 0.327 TRUE 0.667 1.001 8106
## b1 0.931 0.162 0.617 0.928 1.259 FALSE 1.000 1.000 7773
## deviance 74.865 2.006 72.914 74.235 80.188 FALSE 1.000 1.003 2177
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 2 and DIC = 76.872
## DIC is an estimate of expected predictive error (lower is better).Podemos graficar las cadenas y posteriores
plot(jags.sim)
Pregunta:
- Qué tan similares o diferentes son las estimaciones con
JAGScomparadas con nuestro MCMC?
Desafío:
JAGSreporta la devianza (-2 por el logaritmo de likelihood) del modelo. Cómo podemos hacer para obtener la devianza del ajuste de nuestro MCMC?
Stan
Veamos como hacer este mismo análisis en Stan.
El proceso general es parecido a lo que hacemos con JAGS, pero tenemos que definir más detalles en el documento del modelo. Para complicar un poco la cosa, los nombres que se usan para las distribuciones son diferentes a los de JAGS o R . Por suerte, existe un manual muy detallado sobre Stan que lo pueden encontrar aquí.
En el bloque data definimos las cantidades conocidas y qué formato tienen. Luego, definimos en el bloque de parameters a las cantidades no observadas que queremos cuantificar. Finalmente. el bloque model contiene las previas y el modelo de datos.
cat(file = "rp.stan",
"
data {
int<lower=1> n; // total number of observations
int<lower=0> y[n]; // response variable
vector[n] x;
}
parameters {
real b0;
real b1;
}
model{
b0 ~ normal(0,5);
b1 ~ normal(0,1);
y ~ poisson_log(b0 + b1 * x);
}
"
)Para llamar a Stan desde R, usamos el paquete rstan. Para poder correr cadenas en paralelo usamos las opciones rstan_options(auto_write = TRUE) y options(mc.cores = parallel::detectCores()). Luego definimos los datos que vamos a pasarle a Stan y ejecutamos stan definiendo dónde está el modelo, los datos, cuántas iteraciones, el “thining” y cuántas cadenas queremos correr.
library(rstan)## Loading required package: StanHeaders## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:jagsUI':
##
## traceplot## The following object is masked from 'package:coda':
##
## traceplotrstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
rem_dat <- list(n = n,
y = y,
x = x)
stan.sim <- stan(file = 'rp.stan',
data = rem_dat,
iter = 10000,
thin = 1,
chains = 3)Para ver los resultados, podemos hacer un print del objeto de salida de stan
print(stan.sim)## Inference for Stan model: rp.
## 3 chains, each with iter=10000; warmup=5000; thin=1;
## post-warmup draws per chain=5000, total post-warmup draws=15000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## b0 -0.10 0.00 0.23 -0.58 -0.25 -0.09 0.05 0.31 3924 1
## b1 0.93 0.00 0.16 0.62 0.82 0.93 1.04 1.26 3858 1
## lp__ -4.16 0.01 1.01 -6.91 -4.53 -3.85 -3.45 -3.19 4986 1
##
## Samples were drawn using NUTS(diag_e) at Fri Apr 26 12:26:50 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Tenemos varias opciones de gráficos. Veamos algunas:
plot(stan.sim, show_density = TRUE, ci_level = 0.5, fill_color = "purple")## ci_level: 0.5 (50% intervals)## outer_level: 0.95 (95% intervals)
plot(stan.sim, plotfun = "hist", pars = "b0", include = FALSE)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
plot(stan.sim, plotfun = "trace", pars = c("b0", "b1"), inc_warmup = TRUE)
plot(stan.sim, plotfun = "rhat")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Pregunta:
- Cómo compran estas estimaciones con las de
JAGSy las de nuestro MCMC?
brms
También podemos usar el paquete brms de Paul Bürkner que sirve como un atajo para correr modelos en Stan
library(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.8.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:rstan':
##
## loodatos <- data.frame(y, x)
priors <- c(set_prior("normal(0, 1)", class = "b"),
set_prior("normal(0, 5)", class = "Intercept"))
brm.sim <- brm(y ~ x,
data = datos,
family = poisson(),
prior = priors,
iter = 10000,
chains = 3)## Compiling the C++ model## Start samplingPodemos hacer un print para ver los resultados
print(brm.sim)## Family: poisson
## Links: mu = log
## Formula: y ~ x
## Data: datos (Number of observations: 30)
## Samples: 3 chains, each with iter = 10000; warmup = 5000; thin = 1;
## total post-warmup samples = 15000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -0.10 0.23 -0.57 0.32 4415 1.00
## x 0.93 0.16 0.63 1.25 4521 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).y podemos graficar los resultados
plot(brm.sim)
Juan Manuel Morales. April 2019. Last udpared: 2019-04-26