This post is part of a running series on Bayesian MCMC tutorials. For updates, follow @StableMarkets.

## Metropolis Review

Metropolis-Hastings is an MCMC algorithm for drawing samples from a distribution known up to a constant of proportionality, . Very briefly, the algorithm works by starting with some initial draw then running through the following process for times:

- Propose a draw from a given proposal distribution .
- Evaluate posterior up to a proportionality constant under this draw: .
- Evaluate the posterior up to a proportionality constant under the previous draw:
- If , then we
**accept**as a draw from and set - If , then we
**accept**with probability and set or else we**reject**it and retain the previous estimate, .

## Caching

These terms “accept” and “reject” give the impression that we either “keep” or “toss” the results from each iteration. But it’s important to note that while we can toss out if it is rejected, the evaluations and should not be tossed. They should be cached and used for future iterations.

If accepted, then in the next iteration we will need to compare under the new proposal with from this accepted proposal – which we have from the previous iteration. It’s what we used to accept this proposal in the first place.

If rejected, then in the next iteration we will need to compare under a new proposal against from – which we have from the previous iterations.

## Logistic Regression

Caching saves us one evaluation per iteration – cutting complexity by half. To give you a sense of how much computation this saves, consider the case of a logistic regression with binary vector , covariate matrix , and parameter vector .

Where is applied element-wise to yield an vector . For simplicity, let’s say we have an improper, uniform prior. Then,

This is a costly evaluation. Computing requires matrix multiplication – operation. Without caching, we’d have to do this multiplication twice – once under the proposal and once under the current draw.

# Example in R

I wrote a short code implementing both a vanilla Metropolis algorithm (no caching) and a Metropolis with caching for a logistic regression with 50,000 observations and four parameters. I used microbenchmark package to compare run time. The vanilla algorithm’s runtime is about twice as long as the algorithm using caching – as expected.

The relative runtime is actually slightly less than 2 because we still need to do the work of storing the evaluation at each step, but still a huge efficiency gain.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.

Learn more about bidirectional Unicode characters

# full code at | |

# https://github.com/stablemarkets/BayesianTutorials/blob/master/MH_with_caching/ | |

# Author: Arman Oganisian | |

library(microbenchmark) | |

library(LaplacesDemon) | |

library(MASS) | |

source("HelperFunctions.R") # contains mh_vanilla and mh_cache | |

set.seed(10) | |

# hyper-parameters and true values | |

lambda<-c(0,0,0,0) | |

phi<-10 | |

true_beta <- matrix(c(0,2,1,–2),ncol=1) | |

N<-50000 | |

# simulate covariates | |

X1 <- rnorm(N) | |

X2 <- rnorm(N) | |

X3 <- rnorm(N) | |

X <- model.matrix(~ X1 + X2 + X3) | |

# simulate outcome | |

Y <- rbinom(n = N, size = 1, prob = invlogit( X %*% true_beta ) ) | |

## Run benchmark | |

bench<-microbenchmark( | |

# Run Vanialla Metropolis | |

MH_vanilla = mh_vanilla(beta_0 = c(0,0,0,0), # initial value | |

p=4, # number of parameters | |

phi = phi, lambda = lambda, #hyperparameters | |

X = X, Y = Y, # Data | |

#iterations and proposal variance | |

mh_trials=1000, jump_v=.2 ), | |

# Run Metropolis with Cache | |

MH_cache = mh_cache(beta_0 = c(0,0,0,0), | |

phi = phi, lambda = lambda, X = X, Y = Y, | |

mh_trials=1000, jump_v=.2, p=4), | |

times = 10) | |

bench |

There are more efficient ways to make use of the computations of $T(\theta^*)$ and $T(\theta^{(t-1)})$, like Rao-Blackwellisation (Gelfand & Smith, 1990; Casella & Robert, 1996; Douc & Robert, 2011) and prefetching (Banterlé et al. (2015).

I’ve come across Rao-Blackwellization, but not prefetching. Thank you for the reference. By the way, I’ve been following your blog for quite some time and really enjoy your content!