[! Summary] The Metropolis-Hastings algorithm allows us to obtain a sequence of random samples from a probability distribution that is challenging to sample from by generating Markov process which asymptotically reaches a unique stationary distribution. The steps of the algorithm are broadly: (1) set a starting point for the sampling, (2) use a Gaussian random walk to suggest a new point, (3) either accept or reject this new point based on set criteria, and (4) repeat.
1 Description
The Metropolis-Hastings algorithm is a Markov Chain Monte Carlo (MCMC) Methods for obtaining a sequence of random samples from a probability distribution from which direct sampling is difficult. This sequence can be used to approximate the distribution (via a histogram) or to compute an integral (e.g. an expected value).
Like most MCMC algorithms, Metropolis-Hastings is generally used for sampling from multi-dimensional distributions, especially when the number of dimensions is high.
2 Intuition
The Metropolis-Hastings algorithm can draw samples from any probability distribution with probability density \(P(x)\), provided that we know a function \(f(x)\) proportional to the density \(P\) and the values \(f(x)\) can be computed. The proportionality requirement is helpful due to it removing the challenge of computing the normalizing constant.
The algorithm generates a sequence of sample values in such a way that, as more and more sample values are produced, the distribution of values more closely approximates the desired distribution. These sample values are produced iteratively, with the distribution of the next sample being dependent only on the current sample value (thus the sample sequence is a Markov Chain).
Specifically, at each iteration, the algorithm picks a candidate for the next sample based on the current sample value. Then with some probability, the candidate is either accepted (in which case it is used in the next iteration) or it is discarded and the current value is reused.
The probability of acceptance is determined by comparing the values of the function \(f(x)\) of the current and candidate sample values with respect to the desired distribution.
3 Formal Derivation
The purpose of the Metropolis-Hastings algorithm is to generate a collection of states according to a desired distribution \(P(x)\). To accomplish this, the algorithm uses a Markov Process, which asymptotically reaches a unique stationary distribution \(\pi(x)\) such that \(\pi(x)=P(x)\).
A Markov process is uniquely defined by its transition probabilities \(P(x'|x)\), the probability of transitioning from any given state \(x\) to any other given state \(x'\).
The derivation of the algorithm states with the condition of detailed balance (which is a principle of kinetic systems which are decomposed into elementary processes which states that at equilibrium, each elementary process is in equilibrium with its reverse process) \[ P(x'|x)P(x)=P(x|x')P(x'), \] which is rewritten as \[ \frac{P(x'|x)}{P(x|x')}=\frac{P(x')}{P(x)}.\tag{$\star$} \] The approach is to separate the transition in two sub-steps; the proposal and the acceptance-rejection. The proposal distribution \(g(x'|x)\) is the conditional probability of proposing state \(x'\) given \(x\), and the acceptance distribution \(A(x',x)\) is the probability to accept the proposed state \(x'\). The transition probability can be written as the product \[ P(x'|x)=g(x'|x)A(x',x), \] which when substituted into \((\star)\) gives \[ \frac{A(x',x)}{A(x.x')}=\frac{{P(x')g(x|x')}}{P(x)g(x'|x)}. \] Next we choose an acceptance ratio that fulfills the condition above. One common choice is the Metropolis choice: \[ A(x',x)=\min\left( 1, \frac{P(x')}{P(x)} \frac{g(x|x')}{g\left( x'|x \right)} \right). \] For this Metropolis acceptance ratio \(A\), either \(A(x',x)=1\) or \(A(x,x')=1\) and, either way, the condition is satisfied.
4 Algorithm
- Initialize:
- Pick an initial state \(x_0\).
- Set \(t=0\).
- Iterate:
- Generate a random candidate state \(x'\) according to \(g(x'|x_t)\).
- Calculate the acceptance probability \(\frac{A(x',x)}{A(x.x')}=\frac{{P(x')g(x|x')}}{P(x)g(x'|x)}\).
- Accept or Reject:
- Generate a uniform random number \(u\in[0,1]\);
- If \(u\leq A(x',x_t)\), then accept the new state and set \(x_{t+1}=x'\);
- If \(u>A(x',x_{t})\), then reject the new state and copy the old state forward \(x_{{t+1}}=x_{t}\).
- Increment: Set \(t=t+1\).
Provided that specified conditions are met, the empirical distribution of saved states \(x_{0}, \dots, x_{T}\) will approach \(P(x)\). The number of iterations \(T\) required to effectively estimate \(P(x)\) depends on the number of factors, including the relationship between \(P(x)\) and the proposal distribution and the desired accuracy of estimation. For distribution on discrete state spaces, it has to be of the order of the autocorrelation time of the Markov process.
In general it is not clear which distribution \(g(x'|x)\) one should use or the number of iterations necessary for proper estimation; both are free parameters of the method, which must be adjusted to the particular problem at hand.
5 Bayesian Inference
Markov Chain Monte Carlo (MCMC) Methods can be used to draw samples from the posterior distribution of a statistical model. The acceptance probability is given by \[ \mathbb{P}_{{acc}}(\theta_{i}\rightarrow\theta^*)=\min\left( 1, \frac{{L(y|\theta^*)\pi(\theta^*)}}{L(y|\theta_{i})\pi(\theta_{i})}\frac{Q(\theta_{i}|\theta^*)}{Q(\theta^*|\theta_{i})} \right) \] where \(L\) is the likelihood, \(\pi(\theta)\) the prior density and \(Q\) the (conditional) proposal probability.