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
N = 250
K = 3

covariates = replicate(K, rnorm(n = N))
colnames(covariates) = c('X1', 'X2', 'X3')

# create the model matrix with intercept
X = cbind(Intercept = 1, covariates)

# create a normally distributed variable that is a function of the covariates
coefs = c(5, .2, -1.5, .9)
sigma = 2
mu = X %*% coefs
y  = rnorm(N, mu, sigma)

# 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
fit_lm = lm(y ~ ., data = data.frame(X[, -1]))
# 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.

log_posterior <- function(x, y, b, s2) {
  # Args: X is the model matrix; y the target vector; b and s2 the parameters
  # to be estimated

  beta   = b           
  sigma  = sqrt(s2)
  sigma2 = s2
  mu     = X %*% beta
  
  # priors are b0 ~ N(0, sd = 10), sigma2 ~ invGamma(.001, .001)
  priorbvarinv = diag(1/100, 4) 
  prioralpha   = priorbeta = .001
  
  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

update_coef <- function(i, x, y, b, s2) {
  # Args are the same as above but with additional i iterator argument.
  b[i, ] = MASS::mvrnorm(1, mu = b[i-1, ], Sigma = b_var_scale)  # proposal/jumping distribution
  
  # 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
  unidraw = runif(1)
  accept  = unidraw < min(exp(post_diff), 1)  # accept if so
  
  if (accept) b[i,]
  else b[i-1,]
}

# update step for sigma2

update_s2 <- function(i, x, y, b, s2) {
  s2_candidate = rnorm(1, s2[i-1], sd = sigma_scale)
  
  if (s2_candidate < 0) {
    accept = FALSE
  } 
  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])
    
    unidraw = runif(1)
    accept  = unidraw < min(exp(s2_diff), 1)
  }
  
  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. 
nsim   = 5000
warmup = 1000
thin   = 10

b  = matrix(0, nsim, ncol(X))        # initialize beta update matrix
s2 = rep(1, nsim)                    # initialize sigma vector

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.

c_    = 2.4/sqrt(ncol(b))
b_var = vcov(fit_lm) 
b_var_scale  = b_var * c_^2               
sigma_scale  = .9

We can now run and summarize the model with tools from the coda package.

# Run
for (i in 2:nsim) {
  b[i, ] = update_coef(
    i  = i,
    y  = y,
    x  = X,
    b  = b,
    s2 = s2
  )
  
  s2[i] = update_s2(
    i  = i,
    y  = y,
    x  = X,
    b  = b,
    s2 = s2
  )
}

# calculate acceptance rates
b_acc_rate  = mean(diff(b[(warmup+1):nsim,]) != 0)
s2_acc_rate = mean(diff(s2[(warmup+1):nsim]) != 0)         

b_acc_rate
s2_acc_rate

# get final chain
library(coda)

b_mcmc  = as.mcmc(b[seq(warmup + 1,  nsim, by = thin),])
s2_mcmc = as.mcmc(s2[seq(warmup + 1, nsim, by = thin)])

# get summaries
summary(b_mcmc)
summary(s2_mcmc)
b_acc_rate  = mean(diff(b[(warmup+1):nsim,]) != 0)
s2_acc_rate = mean(diff(s2[(warmup+1):nsim]) != 0)
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.

fit_rstan = rstanarm::stan_glm(y ~ ., data = data.frame(X[, -1]))
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