metropolis-hastings.Rmd

# Metropolis Hastings

The following demonstrates a random walk [Metropolis-Hastings
algorithm](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_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.

```{r mh-setup}
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^.

```{r mh-func-posterior}
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.

```{r mh-func-update}
# 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.

```{r mh-initialization1}
# 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.

```{r mh-initialization2}
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 <span class="pack">coda</span> package.

```{r mh-est-code, eval=F}
# 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)
```

```{r mh-est-show, echo=F}
# 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
```


```{r mh-acceptance-rates, echo=1:2}
b_acc_rate  = mean(diff(b[(warmup+1):nsim,]) != 0)
s2_acc_rate = mean(diff(s2[(warmup+1):nsim]) != 0)

kable_df(data.frame(b_acc_rate, s2_acc_rate))
```

Summarize results.

```{r mh-est-summary, echo=F}
# get final chain
library(coda)

b_mcmc = as.mcmc(data.frame(beta    = b[seq(warmup + 1, nsim, by = thin), ], 
                            sigmasq = s2[seq(warmup + 1, nsim, by = thin)]))

# get summaries; compare to lm and stan
# broom.mixed::tidyMCMC(b_mcmc, conf.int = T) %>% 
#   kable_df()
```

The following table is uses <span class="pack" style = "">rstan's</span> <span class="func" style = "">monitor</span> function to produce typical Stan output.

```{r mh-est-monitor, echo=FALSE}
init = rstan::monitor(
  array(
    cbind(b, s2),
    dim = c(5000, 1, 5),
    dimnames = list(NULL, NULL, colnames(b_mcmc))
  ),
  warmup = warmup,
  digits_summary = 5,
  probs  = c(.025, .975),
  se     = FALSE,
  print  = FALSE
) %>% 
  as_tibble(rownames = 'parameter') 

init %>% 
  select(parameter, mean, sd:Rhat, Bulk_ESS, Tail_ESS) %>% 
  kable_df()
```

## Comparison

We can compare to the <span class="func" style = "">lm</span> result or <span class="pack" style = "">rstanarm</span>.

```{r mh-compare, echo=1, results='hide'}
fit_rstan = rstanarm::stan_glm(y ~ ., data = data.frame(X[, -1]))
```

```{r mh-compare-show, echo=FALSE}
tibble(
  parameter = c(names(coef(fit_lm)), 'sigma_sq'),
  fit = init$mean,
  lm  = c(coef(fit_lm), summary(fit_lm)$sigma^2),
  rstanarm  = c(coef(fit_rstan), fit_rstan$stan_summary['sigma','mean']^2)
) %>% 
  kable_df()
```


## Source

Original demo here:

https://m-clark.github.io/bayesian-basics/appendix.html#metropolis-hastings-example