This note serves as an overview of some of the algorithms used in Bayesian statistics, starting with an introduction to Markov chains and Monte Carlo methods; the two fundamental components that make up a Markov Chain Monte Carlo (MCMC) process. First, we talk about the components of an MCMC process: Markov chains and Monte Carlo methods. Then we use these tools understand how the Metropolis-Hastings and Hamiltonian Monte Carlo algorithms explore the posterior distribution.
A random walk is a special case of a Markov chain. It does precisely what you think it does; randomly move through some mathematical space. The canonical example is a particle moving in one dimension. With equal probability the particle moves up or down at time \(t\) relative to it position at time \(t-1\). This behavior can be represented with a Bernoulli probability mass function. Another example would be to let \(\varepsilon ~ \mathcal{N}(0,1)\) and define, \[ X_n = X_{n-1} + \varepsilon \] which results in a continuous one-dimensional random walk.
The code below computes a discrete random walk in two dimensions with equal probabilities being assigned to each movement: up, down, left, or right. It uses the sample()
function but rbinom()
could have been used instead by defining \((1,2,3,4) \equiv (\mbox{up, down, left, right})\).
# RANDOM WALK
loc <- c("up", "down", "left", "right")
dat <- matrix(NA, ncol = 2, nrow = 1e5/2)
dat[1,] <- c(0,0)
for(i in 2:nrow(dat)) {
rand <- sample(loc, 1, prob = rep(0.25,4))
if(rand == "up") {
dat[i,] <- c(0,1) + dat[i-1,]
}
else if(rand == "down") {
dat[i,] <- c(0,-1) + dat[i-1,]
}
else if(rand == "left") {
dat[i,] <- c(1,0) + dat[i-1,]
}
else if(rand == "right") {
dat[i,] <- c(-1,0) + dat[i-1,]
}
else {
dat[i,] <- c(NA,NA)
}
}
Below, the figure on the left plots the first 5,000 steps of the random walk, and the figure on the right plots all 50,000 steps of the walk with the first 5,000 steps highlighted. The darker regions indicate the parts of the process that walked over itself.