# The Metropolis Algorithm

I Know how to apply the Metropolis Algorithm, but I'd be grateful if someone could explain to me the reasoning behind the steps in the algorithm. I've tried in vain looking for the original paper.
Thanks.

-

I'll phrase everything below in very simple terms; if you've already seen some formal probability, measure theory and stochastic methods, feel free to mentally replace the notation by a more formal one.

Let's say you have a distribution of 'configurations' $\phi_k$, each having a probability ('weight') $\pi_k$. Now you're interested in finding the mean value $$\langle O \rangle = \frac{\sum_k O(\phi_k) \pi_k} {\sum_k \pi_k}$$ (replace the sums by integrals if you like). There are a lot of possible reasons for which you'll not be able to calculate both above sums explicitly. In that case, it's advantageous to use a stochastic method, and that's what Monte Carlo does: you're constructing a Markov chain $$\dotsm \rightarrow \phi_{t-1} \rightarrow \phi_t \rightarrow \phi_{t+1} \rightarrow \dotsm$$ where the $t$ is a discrete time variable. Formally, $\phi_t$ is a random variable, but in practice it's the configuration at time $t$.

What you need is the 'transition matrix' $W(\phi_k \rightarrow \phi_n)$ (can be more complicated in general), which describes the probability at each time to go from configuration $\phi_k$ to $\phi_n$. Now one can prove, that if $W$ satisfies a small number of conditions, the above Markov chain actually 'converges', in the sense that as $t \rightarrow \infty,$ the probability to encounter a configuration $\phi_k$ goes to $\pi_k$. The crucial condition for this to happen is $$\sum_k W(\phi_k \rightarrow \phi_n) \pi_k = \sum_k W(\phi_n \rightarrow \phi_k) \pi_n,$$ which expresses that once you're at equilibrium, you'll stay there. One way to satisfy the above condition, is to do it term-by-term, i.e. $$W(\phi_k \rightarrow \phi_n) \pi_k = W(\phi_n \rightarrow \phi_k) \pi_n,$$ called 'detailed balance'. The actual Metropolis formula, $$W(\phi_k \rightarrow \phi_n) = \min(1,\pi_n/\pi_k),$$ is nothing but one of the many solutions to the above condition. Finally, since at $t = \infty$ you're sampling according to the 'right' distribution $\pi_k,$ you can calculate the mean as you go: $$\langle O \rangle \simeq \frac{1}{N} \sum_{t=1}^N O(\phi_t).$$ The error of the above approximation is non-trivial to determine, however.

-
Thanks for your answer. How different is the Metropolis-Hastings Algorithm from the Metropolis Algorithm? –  Nana Apr 13 '11 at 3:25
@Nana: sorry for the delay in my reply! The above supposes that the probability to propose a new configuration is uniform, i.e. if your current configuration is $\phi_k,$ then the probability to consider a new configuration $\phi_{k'}$ is uniform for all $\phi_{k'}.$ If that's not the case, Metropolis-Hastings uses a factor $Q(\phi_k,\phi_{k'})$ to compensate for this. Finally, the nomenclature differs, i.e. some authors always or never use the name Metropolis-Hastings! –  Gerben Apr 16 '11 at 17:11
Thanks –  Nana Apr 17 '11 at 1:40

The idea of the algorithm, is using one distribution (the transition distribution) in order to sample a different one (the original distribution). The assumption is that the original distribution is calculable, but that it is too difficult to sample directly. Ideally, the transition distribution should be "close enough" to the original distribution.

In order to do that, we have to sample in such a way where we don't follow the transition distribution "blindly", but rather use the original distribution to weigh how likely it is to actually sample the next value, rather than the old one.

It can be shown, (by simple summation of the probabilities) that the stationary distribution of the process - if there is one (if I'm not mistaken, in order for there to be one, the transition distribution has to satisfy some conditions that would ensure the sampling is "close enough" to the original distribution) - is indeed the distribution of the original distribution.

A common case which uses this algorithm, is the Gibbs sampling algorithm - on large graphical models. The transition distribution used is a single variable value change.

-
Thanks for the explanation. If you could also explain the mathematics of it, I'd be very grateful. –  Nana Apr 13 '11 at 3:09