Could someone help me to understand the Metropolis-Hasting algorithm for discrete Markov Chains?

Since no answer was given at this question, I am going to make it here.

Metropolis-Hastings Algorithm

Assume the Markov chain is in some state $$X_{n} = i$$. Let $$\textbf{H}$$ be the transition matrix for any irreducible Markov chain on the state space. We generate $$X_{n+1}$$ via the following algorithm:

(a) Choose a proposal state $$j$$ according to the probability distribution in the $$i$$-th row of $$\textbf{H}$$.

(b) Compute the acceptance probability $$\alpha_{ij} = \min\left\{1,\displaystyle\frac{\pi_{j}H_{ji}}{\pi_{i}H_{ij}}\right\}$$

(c) Generate a uniform random number $$U\sim\text{Uniform}(0,1)$$. If $$U<\alpha_{ij}$$, accept the move and set $$X_{n+1} = j$$. Otherwise, reject the move and keep $$X_{n+1} = X_{n}$$.

MY DOUBTS

I know that Markov chains that are aperiodic and irreducible admit a stationary distribution. I also know that any distribution $$\pi$$ satisfying the reversibility condition is an equilibrium distribution.

As far as I have understood, the Metropolis-Hasting algorithm proposes a Markov chain whose stationary distribution is $$\pi$$ (the target distribution we are interested in sampling) and then simulate the Markov chain.

EDIT

However, I am a little bit lost as to the interpretation of the algorithm itself. Could someone interpret it step by step to me?

I know such question is quite naive, but I am not able to understand it properly. Any help is appreciated.

• What do the symbols $\pi_i$ and $\pi_j$ stand for? – Gerry Myerson Oct 7 '19 at 1:15
• I think it is related to the stationary distribution $\pi$ from which we want to sample. – APCorreia Oct 7 '19 at 1:16
• Then perhaps you should edit that into the body of the question. – Gerry Myerson Oct 7 '19 at 1:52
• @user1337 Can you let me know how to improve my answer or, if it is satisfactory, accept it? – Ian Oct 7 '19 at 17:47
• @Ian Can you tell me how do we propose a candidate for the next stage through $\textbf{H}$? How do we utilize $\pi$ in order to do so? – APCorreia Oct 7 '19 at 18:20

The main mathematical content of Metropolis-Hastings follows by examining what its transition probability matrix effectively is, even though this transition probability matrix isn't actually used to run Monte Carlo with it. This matrix has off-diagonal elements $$p_{ij}=\min \left \{ H_{ij},\frac{\pi_j H_{ji}}{\pi_i} \right \}$$ (you just multiply $$H_{ij}$$ by the acceptance probability since you make the move proposal and the move acceptance check independently). Now calculate $$\pi_i p_{ij}=\min \left \{ \pi_i H_{ij},\pi_j H_{ji} \right \}=\pi_j p_{ji}$$. Thus the MH transition probability matrix is reversible with respect to $$\pi$$, so it has $$\pi$$ as its unique stationary distribution.

Regarding the acceptance process, I can tell you that generally speaking, higher acceptance rates contribute to higher convergence rates. But the higher the acceptance rate is, the closer you were to "being done before you started", i.e. the closer you were to having $$H_{ij}=\pi_j$$. So there is a tradeoff going on here. Importance sampling for computing expected values has something very similar: the lowest possible variance is to sample with probability proportional to $$x p(x)$$, but the normalization constant of this distribution is already the expected value you wanted to find.

When $$H$$ is symmetric, the expression for $$\alpha_{ij}$$ simplifies somewhat; the algorithm is sometimes called the Metropolis algorithm in this context. But if $$\pi$$ takes on wildly varying values, then this can cause low acceptance rates and thus slow convergence, so it isn't necessarily an advantage.

One of the advantages of MH is that it does not require us to know $$\pi$$, but only to know $$\frac{\pi_j}{\pi_i}$$. This is often helpful, especially in chains with many states, where we may know $$\pi_i=\frac{q_i}{\sum q_i}$$ but we don't know $$\sum q_i$$. In this situation we can still compute the ratio, because the normalization constants cancel out. A common example is when $$\pi$$ is the Boltzmann distribution for some energy function at some given temperature. In this situation you can define $$q_i=e^{-\frac{E(i)}{k_B T}}$$ and so $$\frac{\pi_j}{\pi_i}=e^{\frac{E(i)-E(j)}{k_B T}}$$.

Regarding your edit, MH as an algorithm says: given a way to simulate from an irreducible* chain $$H$$, a way to compute $$\pi_j/\pi_i$$, and a way to generate Bernoulli random variables, you can sample from the chain with $$p_{ij}=\min \left \{ H_{ij},\frac{\pi_j H_{ji}}{\pi_i} \right \}$$. You do that as follows:

• Propose a candidate for the next state $$j$$ through $$H$$
• Calculate $$\alpha_{ij}=p_{ij}/H_{ij}$$
• Generate a Bernoulli($$\alpha_{ij}$$) random variable and use its result to decide whether to accept the move. This makes it so that the overall probability of a move from $$i$$ to $$j$$ is $$p_{ij}$$.

*Remark: $$H$$ does not need to be aperiodic, but the MH process itself will always be aperiodic.