Metropolis Hastings
The following demonstrates a random walk Metropolis-Hastings algorithm using the data and model from prior sections of the document. I had several texts open while cobbling together this code (noted below), and some oriented towards the social sciences. Some parts of the code reflect information and code examples found therein, and follows Lynch’s code a bit more.
References:
Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed.
Gill, Jeff. 2008. Bayesian Methods : A Social and Behavioral Sciences Approach. Second.
Jackman, Simon. 2009. Bayesian Analysis for the Social Sciences.
Lynch, Scott M. 2007. Introduction to Applied Bayesian Statistics and Estimation for Social Scientists.
Data Setup
Here we create some data based on a standard linear regression.
library(tidyverse)
# set seed for replicability
set.seed(8675309)
# create a N x k matrix of covariates
= 250
N = 3
K
= replicate(K, rnorm(n = N))
covariates colnames(covariates) = c('X1', 'X2', 'X3')
# create the model matrix with intercept
= cbind(Intercept = 1, covariates)
X
# create a normally distributed variable that is a function of the covariates
= c(5, .2, -1.5, .9)
coefs = 2
sigma = X %*% coefs
mu = rnorm(N, mu, sigma)
y
# same as
# y = 5 + .2*X1 - 1.5*X2 + .9*X3 + rnorm(N, mean = 0, sd = 2)
# Run lm for later comparison; but go ahead and examine now if desired
= lm(y ~ ., data = data.frame(X[, -1]))
fit_lm # summary(fit_lm)
Functions
The primary functions that we need to specify regard the posterior distribution, an update step for beta coefficients, and an update step for the variance estimate. We assume a normal distribution for the β coefficients, inverse gamma on σ2.
<- function(x, y, b, s2) {
log_posterior # Args: X is the model matrix; y the target vector; b and s2 the parameters
# to be estimated
= b
beta = sqrt(s2)
sigma = s2
sigma2 = X %*% beta
mu
# priors are b0 ~ N(0, sd = 10), sigma2 ~ invGamma(.001, .001)
= diag(1/100, 4)
priorbvarinv = priorbeta = .001
prioralpha
if (is.nan(sigma) | sigma<=0) { # scale parameter must be positive
return(-Inf)
}
# Note that you will not find the exact same presentation across texts and
# other media for the log posterior in this conjugate setting. In the end
# they are conceptually still (log) prior + (log) likelihood (See commented 'else')
else {
-.5*nrow(X)*log(sigma2) - (.5*(1/sigma2) * (crossprod(y-mu))) +
-.5*ncol(X)*log(sigma2) - (.5*(1/sigma2) * (t(beta) %*% priorbvarinv %*% beta)) +
-(prioralpha + 1)*log(sigma2) + log(sigma2) - priorbeta/sigma2
}# else {
# ll = mvtnorm::dmvnorm(y, mean=mu, sigma=diag(sigma2, length(y)), log=T)
# priorb = mvtnorm::dmvnorm(beta, mean=rep(0, length(beta)), sigma=diag(100, length(beta)), log=T)
# priors2 = dgamma(1/sigma2, prioralpha, priorbeta, log=T)
# logposterior = ll + priorb + priors2
# logposterior
# }
}
Update functions.
# update step for regression coefficients
<- function(i, x, y, b, s2) {
update_coef # Args are the same as above but with additional i iterator argument.
= MASS::mvrnorm(1, mu = b[i-1, ], Sigma = b_var_scale) # proposal/jumping distribution
b[i, ]
# Compare to past- does it increase the posterior probability?
=
post_diff log_posterior(x = x, y = y, b = b[i, ], s2 = s2[i-1]) -
log_posterior(x = x, y = y, b = b[i-1, ], s2 = s2[i-1])
# Acceptance phase
= runif(1)
unidraw = unidraw < min(exp(post_diff), 1) # accept if so
accept
if (accept) b[i,]
else b[i-1,]
}
# update step for sigma2
<- function(i, x, y, b, s2) {
update_s2 = rnorm(1, s2[i-1], sd = sigma_scale)
s2_candidate
if (s2_candidate < 0) {
= FALSE
accept
} else {
=
s2_diff log_posterior(x = x, y = y, b = b[i, ], s2 = s2_candidate) -
log_posterior(x = x, y = y, b = b[i, ], s2 = s2[i - 1])
= runif(1)
unidraw = unidraw < min(exp(s2_diff), 1)
accept
}
ifelse(accept, s2_candidate, s2[i - 1])
}
Estimation
Now we can set things up for the MCMC chain. Aside from the typical MCMC setup and initializing the parameter matrices to hold the draws from the posterior, we also require scale parameters to use for the jumping/proposal distribution. While this code regards only one chain, though a simple loop or any number of other approaches would easily extend it to two or more.
# Setup, starting values etc.
= 5000
nsim = 1000
warmup = 10
thin
= matrix(0, nsim, ncol(X)) # initialize beta update matrix
b = rep(1, nsim) # initialize sigma vector s2
For the following, this c_
term comes from BDA3 12.2 and will produce an acceptance rate of .44 in 1 dimension and declining from there to about .23 in high dimensions. For the sigma_scale, the magic number comes from starting with a value of one and fiddling from there to get around .44.
= 2.4/sqrt(ncol(b))
c_ = vcov(fit_lm)
b_var = b_var * c_^2
b_var_scale = .9 sigma_scale
We can now run and summarize the model with tools from the coda package.
# Run
for (i in 2:nsim) {
= update_coef(
b[i, ] i = i,
y = y,
x = X,
b = b,
s2 = s2
)
= update_s2(
s2[i] i = i,
y = y,
x = X,
b = b,
s2 = s2
)
}
# calculate acceptance rates
= mean(diff(b[(warmup+1):nsim,]) != 0)
b_acc_rate = mean(diff(s2[(warmup+1):nsim]) != 0)
s2_acc_rate
b_acc_rate
s2_acc_rate
# get final chain
library(coda)
= as.mcmc(b[seq(warmup + 1, nsim, by = thin),])
b_mcmc = as.mcmc(s2[seq(warmup + 1, nsim, by = thin)])
s2_mcmc
# get summaries
summary(b_mcmc)
summary(s2_mcmc)
= mean(diff(b[(warmup+1):nsim,]) != 0)
b_acc_rate = mean(diff(s2[(warmup+1):nsim]) != 0) s2_acc_rate
b_acc_rate | s2_acc_rate |
---|---|
0.298 | 0.43 |
Summarize results.
The following table is uses rstan’s monitor function to produce typical Stan output.
parameter | mean | sd | 2.5% | 97.5% | n_eff | Rhat | Bulk_ESS | Tail_ESS |
---|---|---|---|---|---|---|---|---|
beta.1 | 4.900 | 0.135 | 4.638 | 5.162 | 258 | 1.007 | 252 | 235 |
beta.2 | 0.083 | 0.136 | -0.192 | 0.354 | 289 | 1.003 | 294 | 425 |
beta.3 | -1.468 | 0.122 | -1.706 | -1.232 | 304 | 1.001 | 306 | 427 |
beta.4 | 0.831 | 0.120 | 0.598 | 1.061 | 354 | 1.000 | 359 | 518 |
sigmasq | 4.091 | 0.377 | 3.431 | 4.880 | 843 | 1.001 | 832 | 964 |
Comparison
We can compare to the lm result or rstanarm.
= rstanarm::stan_glm(y ~ ., data = data.frame(X[, -1])) fit_rstan
parameter | fit | lm | rstanarm |
---|---|---|---|
(Intercept) | 4.900 | 4.898 | 4.898 |
X1 | 0.083 | 0.084 | 0.088 |
X2 | -1.468 | -1.469 | -1.468 |
X3 | 0.831 | 0.820 | 0.821 |
sigma_sq | 4.091 | 4.084 | 4.110 |
Source
Original demo here:
https://m-clark.github.io/bayesian-basics/appendix.html#metropolis-hastings-example