Markov chains and MCMC: Metropolis–Hastings, Gibbs and Hamiltonian Monte Carlo

Before this lesson: lesson 04 (joint and conditional distributions), lesson 05 (Gaussian conditionals), lesson 10 (posteriors known only up to a constant) and lesson 14 (Monte Carlo estimates, importance sampling and its effective sample size). From linear algebra: left and right eigenvectors (linear algebra 05). Time: about 90 minutes, plus the notebook.

---

Why this matters #

• Bayesian inference. Most posteriors have no closed form (lesson 10). Stan, PyMC and NumPyro sample them with Hamiltonian Monte Carlo (the NUTS variant). This is the standard way to fit hierarchical models, to calibrate the parameters of a physical model with uncertainty (a climate model component, a seismic source), and to get uncertainty for the last layer of a neural network.
• When importance sampling fails. Importance sampling (lesson 14) needs a proposal close to the target. In high dimension the weights collapse onto a few samples. MCMC makes small local moves instead, so it does not need a global proposal.
• Generative models. Energy-based and score-based models sample with Langevin dynamics, a gradient step plus Gaussian noise. This is the one-step version of HMC (generative models 06, generative models 09).
• Markov chains everywhere. PageRank is the stationary distribution of a random surfer. A policy in a Markov decision process turns it into a Markov chain (reinforcement learning 01). An n-gram language model is a Markov chain on words.

1. Predict first #

2. The picture #

A Markov chain moves probability mass. Write the distribution of the state $X_t$ as a row vector $\mu_t$: entry $i$ is $P(X_t = i)$. The transition matrix $P$ has entries $P_{ij} = P(X_{t+1} = j \mid X_t = i)$, so every row sums to 1. One step of the chain moves the mass: $\mu_{t+1} = \mu_t P$. Mass flows out of every state and into every state. The distribution stops changing when, at every state, the inflow equals the outflow. That distribution is the stationary distribution $\pi$. (In this lesson $\pi$ is a distribution; the number $3.14\ldots$ appears only inside formulas such as $\sqrt{2\pi}$.)

Interactive figure: loading…

MCMC reverses the usual question. Usually $P$ is given, and you look for $\pi$. In MCMC, the target $\pi$ is given (up to a constant), and you design $P$ so that $\pi$ is stationary. Then you run the chain and estimate $\mathbb{E}_\pi[f]$ by the time average $\frac1n \sum_t f(X_t)$. The simplest design is random-walk Metropolis: from the current point $x$, propose $x' = x + s z$ with $z \sim \mathcal N(0, 1)$. If the target is higher at $x'$, move. If it is lower, move only with probability $\tilde\pi(x')/\tilde\pi(x)$; otherwise stay at $x$ and record $x$ again. Uphill moves are always taken, downhill moves sometimes. The chain spends more time where the target is high, in exactly the right proportions (section 3 proves this).

Interactive figure: loading…

Gibbs sampling moves along the axes. In two dimensions, Gibbs sampling draws $x_1$ from its conditional distribution given the current $x_2$, then $x_2$ given the new $x_1$. Each move is parallel to an axis, so the path is a staircase. If the target is a narrow diagonal ridge (strong correlation), the steps are short and the chain needs many sweeps to travel along the ridge.

Interactive figure: loading…

HMC slides a puck on a landscape. Turn the target into a surface: the height at $x$ is $U(x) = -\log\tilde\pi(x)$. High probability means low height. Place a frictionless puck at the current point, give it a random kick (a momentum $p$), and let it slide for a while. It speeds up in the valleys and slows down on the slopes, so it follows the shape of the distribution and travels far in one iteration. The computer cannot solve the motion exactly, so it takes small leapfrog steps. A final accept/reject step removes the small error of these steps.

Interactive figure: loading…

3. Definitions and main results #

Finite Markov chains #

A Markov chain on the states $\{1, \dots, K\}$ is a sequence $X_0, X_1, \dots$ with the Markov property:

The next state depends only on the current state. Every row of $P$ sums to 1. By the law of total probability, $\mu_{t+1}(j) = \sum_i \mu_t(i) P_{ij}$, that is $\mu_{t+1} = \mu_t P$, and so $\mu_t = \mu_0 P^t$.

• A stationary distribution is a probability vector $\pi$ with $\pi P = \pi$. So $\pi$ is a left eigenvector of $P$ for the eigenvalue 1 (an eigenvector of $P^\top$). Every finite chain has at least one. The eigenvalue 1 exists because $P\mathbf 1 = \mathbf 1$ (rows sum to 1), and the Perron–Frobenius theorem gives a left eigenvector with nonnegative entries.
• Irreducible: every state can reach every other state in some number of steps. Then $\pi$ is unique and every $\pi_i > 0$. The mean time to return to state $i$ is $1/\pi_i$.
• Aperiodic: the period of state $i$ is the greatest common divisor of the times $t \ge 1$ with $(P^t)_{ii} > 0$. A chain is aperiodic if every period is 1. In an irreducible chain, one state with $P_{ii} > 0$ is enough.

Convergence theorem. If a finite chain is irreducible and aperiodic, then $\mu_0 P^t \to \pi$ for every start $\mu_0$. The error shrinks geometrically, roughly like $|\lambda_2|^t$, where $|\lambda_2| < 1$ is the second-largest absolute value among the eigenvalues of $P$.

Ergodic theorem. If a finite chain is irreducible, then for every start and every function $f$,

This is the law of large numbers for Markov chains, and it is what MCMC uses: time averages along one path estimate expectations under $\pi$. (It holds even for periodic chains, where $\mu_t$ itself does not converge.)

Two states, by hand. Let $P = \begin{pmatrix} 1 - a & a \\ b & 1 - b \end{pmatrix}$ with $a, b \in (0, 1)$.

1. Write $m_t = \mu_t(1)$. Then $m_{t+1} = (1 - a)\,m_t + b\,(1 - m_t) = b + (1 - a - b)\,m_t$.
2. A fixed point satisfies $m = b + (1 - a - b)\,m$, so $m = b/(a + b)$. Hence $\pi = \frac{1}{a + b}(b,\ a)$.
3. Subtract the fixed-point equation from step 1: $m_{t+1} - \pi_1 = (1 - a - b)\,(m_t - \pi_1)$.
4. So the error is multiplied by $\lambda_2 = 1 - a - b$ at every step, from any start. The eigenvalues of $P$ are 1 and $1 - a - b$. $\blacksquare$

Why the conditions matter.

• The flip chain $P = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$ is irreducible with the unique $\pi = (\tfrac12, \tfrac12)$. But it has period 2: started in state 1, $P(X_t = 1)$ alternates between 0 and 1 and never converges. Its eigenvalues are $1$ and $-1$, so $|\lambda_2| = 1$.
• A reducible chain (two groups of states that never communicate) has one stationary distribution for each closed group, and every mixture of them. The long-run behaviour then depends on the start.

Detailed balance #

Definition. A distribution $\pi$ satisfies detailed balance with $P$ if

In words: in the long run, the flow of probability from $i$ to $j$ equals the flow from $j$ to $i$, for every pair. A chain with this property is called reversible: run backwards in time, it has the same law.

Theorem. Detailed balance implies stationarity: $\pi P = \pi$.

Proof. Fix a state $j$.

1. $(\pi P)_j = \sum_i \pi_i P_{ij}$ (definition of the product).
2. By detailed balance, each term $\pi_i P_{ij}$ equals $\pi_j P_{ji}$. So $(\pi P)_j = \sum_i \pi_j P_{ji} = \pi_j \sum_i P_{ji}$.
3. Row $j$ of $P$ sums to 1, so $(\pi P)_j = \pi_j$. $\blacksquare$

Detailed balance is sufficient, not necessary. The Cycle preset in the first figure moves $1 \to 2 \to 3 \to 1$ with probability 0.9 and stays with probability 0.1. Its stationary distribution is uniform, but all flow goes one way around the cycle, so $\pi_1 P_{12} = 0.3 \ne 0 = \pi_2 P_{21}$. MCMC constructions still check detailed balance, because it is local: one pair of states at a time, without a sum over all states.

Metropolis–Hastings #

The target is $\pi(x) = \tilde\pi(x)/Z$. You can evaluate $\tilde\pi$, but you do not know $Z$. For a posterior, $\tilde\pi(\theta) = p(\mathcal D \mid \theta)\,p(\theta)$ and $Z = p(\mathcal D)$, the evidence (lesson 10). A proposal $q(x' \mid x)$ is any distribution you can sample from.

Algorithm (one step). From the current state $x$:

1. propose $x' \sim q(x' \mid x)$;
2. compute the acceptance probability

3. with probability $\alpha$ move to $x'$; otherwise stay at $x$, and record $x$ again.

In code, compare logs: accept if $\log u < \log\tilde\pi(x') - \log\tilde\pi(x) + \log q(x \mid x') - \log q(x' \mid x)$ with $u \sim \text{Uniform}(0, 1)$. This avoids overflow and underflow.

Theorem (Metropolis–Hastings). The Metropolis–Hastings kernel satisfies detailed balance with $\pi$. So $\pi$ is a stationary distribution of the chain.

Proof. We write the proof with densities; for a finite state space replace integrals by sums.

1. For $x' \ne x$, the chain moves from $x$ to $x'$ with density $K(x, x') = q(x' \mid x)\,\alpha(x, x')$. On a rejection it stays at $x$; this adds a point mass at $x' = x$.
2. Multiply by $\pi(x)$ and move the factor inside the minimum:

The right side does not change if we swap $x$ and $x'$. So $\pi(x)K(x, x') = \pi(x')K(x', x)$ for $x' \ne x$. 3. The point mass at $x' = x$ satisfies detailed balance trivially: both sides are the same term. 4. Integrate detailed balance over $x$: $\int \pi(x)K(x, x')\,dx = \pi(x') \int K(x', x)\,dx = \pi(x')$, because the kernel from $x'$, including its rejection mass, has total probability 1. This is stationarity. $\blacksquare$

Only ratios appear. $\tilde\pi(x')/\tilde\pi(x) = \pi(x')/\pi(x)$: the constant $Z$ cancels. This is why MCMC works for posteriors. Stationarity alone does not give convergence: the chain must also be irreducible and aperiodic. With a Gaussian proposal on $\mathbb R^d$ and a positive target, both hold (rejections give aperiodicity).

• Random-walk Metropolis. $x' = x + s z$ with $z \sim \mathcal N(0, I)$. This proposal is symmetric, $q(x' \mid x) = q(x \mid x')$, so $\alpha = \min(1, \tilde\pi(x')/\tilde\pi(x))$. A small step $s$: almost every move is accepted, but the moves are tiny. A large $s$: almost every move is rejected. For a 1-D $\mathcal N(0, 1)$ target the acceptance rate is $\frac{2}{\pi}\arctan(2/s)$: 0.844, 0.442 and 0.205 for $s = 0.5, 2.4, 6$. In $d$ dimensions the best step scales like $2.38/\sqrt d$, with an acceptance rate near 0.234 (Roberts, Gelman and Gilks, 1997), and the chain needs about $O(d)$ steps per independent sample.
• The Hastings correction $q(x \mid x')/q(x' \mid x)$ is needed when the proposal is not symmetric. Example: for a positive parameter, the multiplicative walk $x' = x\,e^{s z}$ has the correction $x'/x$ (check yourself, question 4).
• Independence sampler. The proposal $q(x')$ does not depend on $x$. Then $\alpha = \min\big(1, w(x')/w(x)\big)$ with the importance weight $w = \tilde\pi/q$ of lesson 14. As in importance sampling, $q$ needs heavier tails than $\pi$.

Gibbs sampling #

Let $x = (x_1, \dots, x_d)$. The full conditional of coordinate $i$ is $\pi(x_i \mid x_{-i})$, where $x_{-i}$ is the vector of all other coordinates. Gibbs sampling updates one coordinate at a time with a draw from its full conditional. A sweep updates every coordinate once.

Theorem. One Gibbs update is a Metropolis–Hastings step with acceptance probability 1.

Proof. The proposal changes only coordinate $i$: $q(x' \mid x) = \pi(x'_i \mid x_{-i})$, with $x'_{-i} = x_{-i}$.

1. Write $\pi(x) = \pi(x_{-i})\,\pi(x_i \mid x_{-i})$. Then $\pi(x)\,q(x' \mid x) = \pi(x_{-i})\,\pi(x_i \mid x_{-i})\,\pi(x'_i \mid x_{-i})$.
2. In the same way, $\pi(x')\,q(x \mid x') = \pi(x_{-i})\,\pi(x'_i \mid x_{-i})\,\pi(x_i \mid x_{-i})$, using $x'_{-i} = x_{-i}$.
3. The two products are equal, so the MH ratio is 1 and $\alpha = 1$. $\blacksquare$

Each update leaves $\pi$ invariant, so a full sweep (one update after another) does too. A sweep in fixed order is in general not reversible, but it is still a valid sampler. Gibbs sampling is popular when the full conditionals have closed forms, as in conjugate hierarchical models, mixture models and topic models.

Why Gibbs is slow for correlated targets. Take $(X_1, X_2)$ standard Gaussian with correlation $\rho$. By lesson 05, $X_1 \mid X_2 = x_2 \sim \mathcal N(\rho x_2,\ 1 - \rho^2)$, and the same with the roles swapped.

1. One sweep: $x_1^{(t+1)} = \rho\,x_2^{(t)} + \sqrt{1 - \rho^2}\,e_1$, then $x_2^{(t+1)} = \rho\,x_1^{(t+1)} + \sqrt{1 - \rho^2}\,e_2$, with independent standard normal noise $e_1, e_2$.
2. Substitute the second line of the previous sweep into the first: $x_1^{(t+1)} = \rho^2 x_1^{(t)} + \text{noise}$. So the chain of $x_1$ is an AR(1) process with coefficient $\varphi = \rho^2$.
3. Its autocorrelation at lag $t$ is $\varphi^t = \rho^{2t}$. The next section turns this into an effective sample size: $n\,(1 - \rho^2)/(1 + \rho^2)$.

The picture: Gibbs moves parallel to the axes, and in a narrow diagonal ridge each move has length of order $\sqrt{1 - \rho^2}$. The fixes are to update correlated coordinates together (blocking) or to change the coordinates (reparametrisation), for example rotate them so that they are uncorrelated.

Hamiltonian Monte Carlo #

Momentum. Add an auxiliary variable $p \in \mathbb R^d$, the momentum, with $p \sim \mathcal N(0, I)$, independent of $x$. The joint density is $\propto e^{-H(x, p)}$ with the Hamiltonian (total energy)

$U$ is the potential energy and $\tfrac12\|p\|^2$ the kinetic energy. The $x$-marginal of $e^{-H}$ is $\pi$: sample $(x, p)$ and drop $p$. (Implementations often use a mass matrix $M$: $p \sim \mathcal N(0, M)$ and kinetic energy $\tfrac12 p^\top M^{-1} p$. Here $M = I$.)

Exact dynamics. Hamilton's equations are $\dot x = p$ and $\dot p = -\nabla U(x)$. Their flow has three properties:

1. it conserves $H$ (the derivative of $H$ along the flow is $\nabla U \cdot p + p \cdot (-\nabla U) = 0$);
2. it preserves volume in $(x, p)$ space (Liouville's theorem);
3. it is reversible: flip the momentum at the end, run for the same time, and you return to the start.

Together these mean the flow leaves $e^{-H}$ invariant. Because $H$ is conserved, one trajectory stays on one energy level. A fresh momentum at every iteration (a Gibbs step on $p$) moves the chain between levels.

Leapfrog. The exact flow is not available, so we discretise it. With step size $\varepsilon$, one leapfrog step is

Consecutive half steps in $p$ merge, so $L$ steps cost $L + 1$ gradient evaluations.

• Volume preserving. Each of the three updates is a shear: it changes one block of variables by an amount that depends only on the other block. Its Jacobian matrix is block triangular with identity blocks on the diagonal, so its determinant is 1. The determinant of the composition is the product, also 1.
• Reversible. Start from the end point with the momentum flipped, $(x^*, -p^*)$, and run $L$ steps: you return to $(x, -p)$. The three updates undo each other in reverse order.
• Accurate. The energy error over a fixed time is $O(\varepsilon^2)$, and it does not grow over long trajectories when $\varepsilon$ is below a stability limit. Halving $\varepsilon$ cuts the error by about 4. Above the limit (for a Gaussian, $\varepsilon > 2\sigma_{\min}$, twice the smallest standard deviation) the path explodes.

HMC (one iteration). Draw $p \sim \mathcal N(0, I)$. Run $L$ leapfrog steps from $(x, p)$ to $(x^*, p^*)$. Accept $x^*$ with probability

otherwise keep $x$. This is a Metropolis–Hastings step on $(x, p)$ with the proposal $(x^*, -p^*)$. The momentum flip makes the proposal map its own inverse, so it is symmetric. Volume preservation means no Jacobian factor appears. Without the correction, the $O(\varepsilon^2)$ energy error would bias the samples.

Why HMC scales. Random-walk Metropolis needs about $O(d)$ steps per independent sample in $d$ dimensions. HMC uses the gradient to make long, directed moves. For i.i.d. product targets the cost grows like $O(d^{1/4})$ leapfrog steps per sample, with the best acceptance rate near 0.65 (Beskos et al., 2013). The tuning parameters are $\varepsilon$, $L$ and $M$. NUTS (the No-U-Turn Sampler, the default in Stan) chooses the trajectory length automatically, and tunes $\varepsilon$ and $M$ during warm-up. A divergence is a trajectory whose energy error explodes. It warns that the chain cannot explore a region of high curvature, and the results may be biased.

Langevin dynamics, MALA and Brownian motion

HMC with a single leapfrog step ($L = 1$) and a fresh momentum gives the update $x' = x + \tfrac{\varepsilon^2}{2}\nabla\log\tilde\pi(x) + \varepsilon z$ with $z \sim \mathcal N(0, I)$. Write $\eta = \varepsilon^2$: $x' = x + \tfrac{\eta}{2}\nabla\log\pi(x) + \sqrt{\eta}\,z$. This is the Euler step of the Langevin diffusion $dX_t = \tfrac12\nabla\log\pi(X_t)\,dt + dB_t$, where $B_t$ is Brownian motion: the limit of a random walk whose steps scale like $\sqrt{\Delta t}$ (lesson 06, the deeper block on random walks). The Langevin diffusion has $\pi$ as its stationary distribution. With a Metropolis–Hastings correction the update is MALA (Metropolis-adjusted Langevin algorithm). Without the correction, and with mini-batch gradients, it is SGLD (stochastic gradient Langevin dynamics). Score-based generative models run the same update with a learned score $\nabla\log p$ in place of the true one. SDEs and their stationary distributions are in ML for science 12; score-based diffusion samplers in generative models 09.

Diagnostics: effective sample size and split-$\hat R$ #

MCMC draws are correlated, and the chain needs time to forget its start. Two numbers tell you whether to trust the output.

Effective sample size (ESS). Let the chain be stationary with variance $\sigma^2$ and autocorrelation $\rho_t$ at lag $t$ (the correlation between $x_s$ and $x_{s+t}$). Under mild conditions the Markov chain central limit theorem gives

$\tau$ is the integrated autocorrelation time. The ESS is $n/\tau$: your $n$ correlated draws give the same error as $n/\tau$ independent ones. The Monte Carlo standard error of the mean is $\sigma/\sqrt{\text{ESS}}$.

For an AR(1) chain with coefficient $\varphi$: $\rho_t = \varphi^t$, so $\tau = 1 + 2\varphi/(1 - \varphi) = (1 + \varphi)/(1 - \varphi)$ and $\text{ESS}/n = (1 - \varphi)/(1 + \varphi)$. For $\varphi = 0.9$, $\tau = 19$. For $\varphi = -0.5$, $\text{ESS} = 3n$: negatively correlated draws can beat independent ones.

Estimating $\tau$. The sample autocorrelations $\hat\rho_t$ at large lags are mostly noise, so the full sum is useless. Geyer's initial positive sequence adds pairs $\Gamma_k = \hat\rho_{2k} + \hat\rho_{2k+1}$ for $k = 0, 1, \dots$ while they are positive, and stops at the first $\Gamma_k \le 0$. Then $\hat\tau = -1 + 2\sum_k \Gamma_k$. (The $-1$ removes the double-counted $\hat\rho_0 = 1$.)

Split-$\hat R$. Run $m$ chains of length $n$ from dispersed starts. Split each chain into halves, to get $2m$ half-chains of length $h = n/2$.

1. $W$ = the mean of the $2m$ within-chain sample variances.
2. $B = h \times$ the sample variance of the $2m$ half-chain means.
3. $\widehat{\operatorname{var}}^+ = \frac{h - 1}{h}\,W + \frac{B}{h}$ and $\hat R = \sqrt{\widehat{\operatorname{var}}^+ / W}$.

If all chains explore the same distribution, $B/h$ is small and $\hat R \approx 1$. Chains stuck in different regions make $B$ large, so $\hat R > 1$. The split also catches a chain that drifts: its two halves disagree, even when all chains drift the same way. Vehtari et al. (2021) recommend at least 4 chains, $\hat R < 1.01$, and an ESS above 400 (in total over the chains) before you trust the estimates.

Interactive figure: loading…

What diagnostics cannot see. If all chains are stuck in the same mode of a multimodal target, they agree and $\hat R \approx 1$. A region that no chain reaches (for example the narrow neck of a funnel-shaped posterior) does not show up in the ESS. So: run several chains from dispersed starts, check divergences for HMC, and compare with what you know about the model.

What to report. The sampler and its settings, the number of chains, the warm-up length and the number of kept draws, $\hat R$ and ESS (bulk and tail) for every reported quantity, the number of divergences, and the Monte Carlo standard error of each estimate.

Rank-normalised $\hat R$, bulk and tail ESS

Classic $\hat R$ uses means and variances, so it fails for heavy-tailed targets (infinite variance) and misses chains that differ only in their spread. Vehtari et al. (2021) first replace all draws by their normal scores: pool the draws of all chains, rank them, and map rank $r$ to $\Phi^{-1}\big((r - 3/8)/(S + 1/4)\big)$, with $S$ the total number of draws. Split-$\hat R$ on these scores is the bulk-$\hat R$. The same on the folded draws $|x - \text{median}|$ checks the spread (tails). The reported $\hat R$ is the maximum of the two. Bulk ESS is the ESS of the normal scores; tail ESS is the smaller ESS of the indicators $\mathbb 1\{x \le q_{0.05}\}$ and $\mathbb 1\{x \le q_{0.95}\}$. ArviZ, Stan and PyMC report these by default.

Answers to section 1.

1. Two thirds. With $a = 0.2$ (dry to wet) and $b = 0.4$ (wet to dry), $\pi = (b, a)/(a + b) = (2/3, 1/3)$. The chain is irreducible and aperiodic, so day 0 does not matter in the long run. Its influence decays like $0.4^t$.
2. About 100. The chain of $x_1$ is AR(1) with $\varphi = \rho^2 = 0.9801$. So $\text{ESS}/n = (1 - 0.9801)/(1 + 0.9801) = 0.0199/1.9801 \approx 1/99.5$.

4. Worked example #

(a) The weather chain. $P = \begin{pmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{pmatrix}$, state 1 = dry, state 2 = wet.

• Stationary distribution: $\pi = (2/3,\ 1/3)$. Check the first entry of $\pi P$: $\tfrac23 \cdot 0.8 + \tfrac13 \cdot 0.4 = \tfrac{1.6 + 0.4}{3} = \tfrac23$. ✓
• Detailed balance: $\pi_1 P_{12} = \tfrac23 \cdot 0.2 = \tfrac{0.4}{3}$ and $\pi_2 P_{21} = \tfrac13 \cdot 0.4 = \tfrac{0.4}{3}$. ✓ (Every two-state chain with $\pi > 0$ is reversible: there is only one pair, and stationarity forces the two flows to be equal.)
• Speed: $\lambda_2 = 1 - 0.2 - 0.4 = 0.4$. Starting dry, $P(\text{dry on day } t) = \tfrac23 + \tfrac13 \cdot 0.4^t$. This gives $0.8$, $0.72$ and $0.6701$ for $t = 1, 2, 5$.
• Mean return times $1/\pi_i$: 1.5 days back to a dry day, 3 days back to a wet day.

(b) The Metropolis kernel as a matrix. Target $\tilde\pi = (1, 2, 4, 2, 1)$ on states $1, \dots, 5$, so $\pi = (1, 2, 4, 2, 1)/10$. Proposal: left or right neighbour, probability $\tfrac12$ each; a move off the ends is rejected.

• From state 3 ($\tilde\pi = 4$): to state 2, $\alpha = \min(1, 2/4) = \tfrac12$, so $P_{32} = \tfrac12 \cdot \tfrac12 = 0.25$. In the same way $P_{34} = 0.25$, and $P_{33} = 1 - 0.5 = 0.5$.
• From state 2 ($\tilde\pi = 2$): to state 3, $\alpha = \min(1, 4/2) = 1$, so $P_{23} = 0.5$; to state 1, $\alpha = \tfrac12$, so $P_{21} = 0.25$; and $P_{22} = 0.25$.
• Detailed balance for the pair (2, 3): $\pi_2 P_{23} = 0.2 \cdot 0.5 = 0.1$ and $\pi_3 P_{32} = 0.4 \cdot 0.25 = 0.1$. ✓
• The eigenvalues of this $5 \times 5$ kernel are $1, 0.75, 0.5, 0.25, 0$ (section 5 checks them), so the error shrinks like $0.75^t$.

(c) Three steps of random-walk Metropolis by hand. Target $\tilde\pi(x) = e^{-x^2/2}$, start $x_0 = 0$, step $s = 1$, noise $z = (1.2, -0.4, 0.9)$ and uniforms $u = (0.3, 0.9, 0.5)$. Accept if $u < \alpha$.

step	$x$	$x' = x + z$	$\log\frac{\tilde\pi(x')}{\tilde\pi(x)} = -\frac{x'^2 - x^2}{2}$	$\alpha$	$u$	result
1	0	1.2	$-0.72$	0.487	0.3	accept: 1.2
2	1.2	0.8	$+0.40$	1	0.9	accept: 0.8
3	0.8	1.7	$-1.125$	0.325	0.5	reject: 0.8

The chain is $(1.2,\ 0.8,\ 0.8)$. The rejected step records the current state again. Dropping repeated states would bias the chain toward regions where proposals are often accepted.

(d) One HMC iteration by hand. Target $\mathcal N(0, 1)$, so $U(x) = x^2/2$ and $\nabla U(x) = x$. Start at $x = 0$ with the momentum draw $p = 1$, step $\varepsilon = 0.5$, one leapfrog step ($L = 1$).

1. Energy at the start: $H = \tfrac12 \cdot 0^2 + \tfrac12 \cdot 1^2 = 0.5$.
2. Half step in momentum: $p = 1 - 0.25 \cdot 0 = 1$.
3. Full step in position: $x = 0 + 0.5 \cdot 1 = 0.5$.
4. Half step in momentum: $p = 1 - 0.25 \cdot 0.5 = 0.875$.
5. Energy at the end: $H = \tfrac12 \cdot 0.25 + \tfrac12 \cdot 0.765625 = 0.5078$. The energy error is $\Delta H = 0.0078$.
6. Accept with probability $\min(1, e^{-0.0078}) = 0.992$.

The exact flow would conserve $H$ and give $\alpha = 1$. The leapfrog error is small, so almost every proposal is accepted, and the accept/reject step removes the remaining bias.

(e) Split-$\hat R$ catches a drift. One chain with draws $(0, 1, 2, 3)$. The halves are $(0, 1)$ and $(2, 3)$, so $h = 2$.

1. Within-half variances: $0.5$ and $0.5$, so $W = 0.5$.
2. Half means $0.5$ and $2.5$; their sample variance is $2$, so $B = 2 \cdot 2 = 4$.
3. $\widehat{\operatorname{var}}^+ = \tfrac12 \cdot 0.5 + 4/2 = 2.25$ and $\hat R = \sqrt{2.25/0.5} = \sqrt{4.5} \approx 2.12$.

Without the split there is only one chain, no between-chain variance and no warning.

5. In code #

NumPy checks the hand computations. The notebook for this lesson has these cells and the exercises.

import numpy as np
np.set_printoptions(precision=4, suppress=True)

# (a) weather chain: stationary distribution by a linear solve, and convergence
P = np.array([[0.8, 0.2], [0.4, 0.6]])
A = np.vstack([(P.T - np.eye(2))[:-1], np.ones(2)])   # pi (P - I) = 0, with one equation replaced by sum = 1
pi = np.linalg.solve(A, [0.0, 1.0])
print(pi, np.linalg.matrix_power(P, 5)[0])             # [0.6667 0.3333] [0.6701 0.3299]
print(np.sort(np.linalg.eigvals(P).real))              # [0.4 1. ]

# (b) the Metropolis kernel on 5 states is a matrix: check detailed balance
pt = np.array([1, 2, 4, 2, 1.0])                       # unnormalised target
K = len(pt)
Pmh = np.zeros((K, K))
for i in range(K):
    for j in (i - 1, i + 1):
        if 0 <= j < K:                                 # a move off the ends is rejected
            Pmh[i, j] = 0.5 * min(1.0, pt[j] / pt[i])
    Pmh[i, i] = 1.0 - Pmh[i].sum()                     # rejection mass stays at i
pi5 = pt / pt.sum()
flows = pi5[:, None] * Pmh                             # flows[i, j] = pi_i P_ij
print(np.allclose(flows, flows.T), np.allclose(pi5 @ Pmh, pi5))   # True True
print(np.sort(np.abs(np.linalg.eigvals(Pmh)))[::-1])   # [1.   0.75 0.5  0.25 0.  ]

# (c) three random-walk Metropolis steps with fixed noise
log_t = lambda x: -0.5 * x**2
x, chain = 0.0, []
for z, u in zip([1.2, -0.4, 0.9], [0.3, 0.9, 0.5]):
    xp = x + z
    if np.log(u) < log_t(xp) - log_t(x):               # log acceptance test
        x = xp
    chain.append(x)
print(np.round(chain, 2))                              # [1.2 0.8 0.8]

# (d) one leapfrog step for U(x) = x^2/2, then the MH correction
x, p, eps = 0.0, 1.0, 0.5
H = lambda x, p: 0.5 * x**2 + 0.5 * p**2
H0 = H(x, p)
p = p - eps / 2 * x; x = x + eps * p; p = p - eps / 2 * x
print(x, p, H(x, p) - H0, np.exp(H0 - H(x, p)))        # 0.5 0.875 0.0078125 0.9922...

# (e) split-R-hat of the drifting chain (0, 1, 2, 3)
halves = np.array([[0.0, 1.0], [2.0, 3.0]])
h = 2
W = halves.var(axis=1, ddof=1).mean()
B = h * halves.mean(axis=1).var(ddof=1)
print(W, B, np.sqrt(((h - 1) / h * W + B / h) / W))  # 0.5 4.0 2.1213...

A simulation checks the acceptance-rate formula of random-walk Metropolis and the AR(1) behaviour of Gibbs ($\rho = 0.99$):

import numpy as np

def rwm(log_p, x0, n, s, rng):
    x, lp, out, acc = x0, log_p(x0), np.empty(n), 0
    for t in range(n):
        xp = x + s * rng.standard_normal()
        lpp = log_p(xp)
        if np.log(rng.random()) < lpp - lp:
            x, lp, acc = xp, lpp, acc + 1
        out[t] = x                                     # a rejection records x again
    return out, acc / n

rng = np.random.default_rng(0)
for s in (0.5, 2.4, 6.0):
    chain, acc = rwm(lambda x: -0.5 * x**2, 0.0, 50_000, s, rng)
    print(s, round(acc, 3), round(2 / np.pi * np.arctan(2 / s), 3))
# about: 0.5 0.842 0.844 | 2.4 0.445 0.442 | 6.0 0.203 0.205  (Monte Carlo error about ±0.004)

rho, n = 0.99, 200_000
x1, x2, xs = 0.0, 0.0, np.empty(n)
sd = np.sqrt(1 - rho**2)
z = rng.standard_normal((n, 2))
for t in range(n):
    x1 = rho * x2 + sd * z[t, 0]                       # x1 | x2 ~ N(rho x2, 1 - rho^2)
    x2 = rho * x1 + sd * z[t, 1]                       # x2 | x1, with the new x1
    xs[t] = x1
print(round(np.corrcoef(xs[:-1], xs[1:])[0, 1], 3), rho**2)   # about 0.979, and 0.9801

6. Where it shows up #

• Probabilistic programming. Stan, PyMC and NumPyro run NUTS by default: 4 chains, a warm-up that tunes $\varepsilon$ and the mass matrix, then $\hat R$, bulk and tail ESS and the divergence count in the summary table. Typical uses: hierarchical models for A/B tests over many groups, epidemiological models, calibration of physical model parameters (climate model components, a seismic source inversion), pharmacokinetics.
• Gibbs sampling with conjugate conditionals. Latent Dirichlet allocation (topic models), Bayesian mixture models, Bayesian linear regression with a conjugate prior, and data augmentation for probit regression all have closed-form full conditionals.
• Langevin dynamics in deep learning. Energy-based models and score-based diffusion models generate samples with (annealed) Langevin steps; SGLD gives approximate posterior samples of neural network weights from mini-batch gradients (generative models 06).
• Markov chains as models. PageRank is the stationary distribution of a random surfer with random restarts (the restarts make the chain irreducible and aperiodic). A fixed policy turns an MDP into a Markov chain, and its stationary distribution weights the states in policy-gradient theorems (reinforcement learning 01). Hidden Markov models (speech, genomics) put a Markov chain under noisy observations.
• Statistical physics. The original Metropolis algorithm (1953) sampled the Boltzmann distribution $e^{-E/kT}$ of particle systems. The Ising model and molecular simulations still use Metropolis, Gibbs (heat-bath) and Hamiltonian (molecular dynamics) moves.

7. Common mistakes #

• Computing $\pi$ as a high power of $P$. It fails for periodic chains and is slow when $|\lambda_2|$ is close to 1. Solve $\pi(P - I) = 0$ with $\sum_i \pi_i = 1$.
• Using a right eigenvector instead of a left one. $P\mathbf 1 = \mathbf 1$ always holds, so the right eigenvector for eigenvalue 1 is the constant vector. $\pi$ is the eigenvector of $P^\top$.
• "Detailed balance is necessary for stationarity." It is only sufficient (the cycle chain).
• Forgetting the Hastings correction for an asymmetric proposal, for example a random walk on $\log x$ for a positive parameter (the factor $x'/x$).
• Not recording the current state after a rejection, or re-proposing until a move is accepted. The repeated values are part of the chain.
• Tuning for a very high acceptance rate. Near 100% acceptance usually means tiny steps and slow mixing.
• Adapting the step size forever. If the proposal keeps changing with the chain's history, the chain may no longer have $\pi$ as its stationary distribution. Adapt during warm-up only, then fix the settings.
• Thinning to "fix" autocorrelation. Keeping every $k$-th draw throws information away. It only saves memory.
• Trusting one chain, or the old threshold $\hat R < 1.1$. Run at least 4 chains from dispersed starts and use $\hat R < 1.01$.
• Ignoring divergences in HMC. Even a few divergences mean that part of the posterior was not explored.

8. Check yourself #

1. Find $\pi$ for $P = \begin{pmatrix} 0.5 & 0.5 & 0 \\ 0.25 & 0.5 & 0.25 \\ 0 & 0.5 & 0.5 \end{pmatrix}$ without solving a $3 \times 3$ system.
   Answer

   The chain moves only to neighbours, so try detailed balance: $\pi_1 \cdot 0.5 = \pi_2 \cdot 0.25$ gives $\pi_2 = 2\pi_1$, and $\pi_2 \cdot 0.25 = \pi_3 \cdot 0.5$ gives $\pi_3 = \pi_1$. Normalise: $\pi = (1/4,\ 1/2,\ 1/4)$. Detailed balance implies stationarity, so this is $\pi$.
2. A chain has $P = \begin{pmatrix} 0.9 & 0.1 \\ 0.5 & 0.5 \end{pmatrix}$. Give $\pi$, $\lambda_2$, and the number of steps until the distance to $\pi$ is below 1% of its starting value.
   Answer

   $a = 0.1$, $b = 0.5$, so $\pi = (0.5, 0.1)/0.6 = (5/6,\ 1/6)$ and $\lambda_2 = 1 - 0.6 = 0.4$. The error is multiplied by $0.4$ each step: $0.4^t < 0.01$ for $t > \ln 100/\ln 2.5 \approx 5.03$, so 6 steps.
3. A Metropolis–Hastings sampler uses the symmetric proposal $x' = x + s z$ on a target with $\tilde\pi(x) = 0$ for $x < 0$. The chain is at $x = 0.1$ and proposes $x' = -0.3$. What happens, and what goes wrong if you "reflect" the proposal to $+0.3$ instead?
   Answer

   $\tilde\pi(x') = 0$, so $\alpha = 0$: the move is rejected and the chain records $0.1$ again. Reflecting $x' \mapsto |x'|$ is also a valid symmetric proposal (the density of $|x + sz|$ at $y$ is symmetric in $x$ and $y$ for $x, y \ge 0$), so it is not wrong. What is wrong is to reject and then re-propose until a positive value appears: that changes the kernel near 0 and biases the samples.
4. For a positive parameter you propose $x' = x\,e^{s z}$ with $z \sim \mathcal N(0, 1)$. What is the Hastings correction?
   Answer

   Given $x$, $x'$ is log-normal: $q(x' \mid x) = \frac{1}{x' s\sqrt{2\pi}}\exp\!\big(-\frac{(\log x' - \log x)^2}{2s^2}\big)$. The exponent is symmetric in $x$ and $x'$, so $q(x \mid x')/q(x' \mid x) = x'/x$, and $\alpha = \min\big(1, \frac{\tilde\pi(x')\,x'}{\tilde\pi(x)\,x}\big)$. Equivalently, run a symmetric random walk on $y = \log x$ with the target $\tilde\pi(e^y)\,e^y$ (the Jacobian of lesson 05).
5. Why does the HMC acceptance probability contain no Jacobian term, and why is a fresh momentum drawn at every iteration?
   Answer

   The leapfrog map is volume preserving (a composition of shears with determinant 1), so the density of the proposed point changes only through $e^{-H}$. With the momentum flip the map is its own inverse, so the proposal is symmetric. A fresh momentum is needed because $H$ is (nearly) conserved along a trajectory: without it the chain stays on one energy level and cannot explore the whole target.
6. Four chains give $\hat R = 1.00$ and an ESS of 3000 for every parameter. Name one way the result can still be wrong.
   Answer

   All four chains can sit in the same mode of a multimodal posterior (they agree with each other), or all can miss a hard region such as the neck of a funnel. $\hat R$ and ESS only see the regions the chains visited. For HMC, check divergences; run more chains from more dispersed starts.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What equation defines a stationary distribution $\pi$ of a finite Markov chain with transition matrix $P$?

   Show answer

2. Which two conditions make a finite chain converge to a unique $\pi$ from every start, and how fast?

   Show answer

3. State detailed balance and what it implies.

   Show answer

4. Compute $\pi$ for $P = \begin{pmatrix} 0.8 & 0.2 \\ 0.4 & 0.6 \end{pmatrix}$.

   Show answer

5. Write the Metropolis–Hastings acceptance probability.

   Show answer

6. Why does Metropolis–Hastings not need the normalising constant $Z$?

   Show answer

7. Why is a Gibbs update accepted with probability 1?

   Show answer

8. Why is Gibbs slow for a bivariate Gaussian with $\rho = 0.99$?

   Show answer

9. Name three properties of the leapfrog integrator that HMC needs.

   Show answer

10. How is the effective sample size computed from autocorrelations?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you sample from an unnormalised posterior?" — Metropolis–Hastings: the acceptance probability uses only $\tilde\pi(x')/\tilde\pi(x)$, so $Z$ cancels. Tune the step size, run several chains, check $\hat R$ and ESS. In practice, use NUTS in Stan, PyMC or NumPyro.
• "Prove that Metropolis–Hastings samples the target." — Show detailed balance: $\pi(x)q(x' \mid x)\alpha(x, x') = \min(\pi(x)q(x' \mid x), \pi(x')q(x \mid x'))$ is symmetric; summing gives stationarity; irreducibility and aperiodicity give convergence.
• "Gibbs, random-walk Metropolis or HMC?" — Gibbs when the full conditionals are easy (conjugate models), but it is slow for correlated coordinates. HMC for smooth, continuous, high-dimensional targets with gradients. Random-walk Metropolis for low dimension or when no gradient is available. Discrete parameters: Gibbs or marginalise them out.
• "Importance sampling or MCMC for a posterior?" — Importance sampling in low dimension with a good heavy-tailed proposal: independent draws and an estimate of the evidence $Z$. MCMC in high dimension, where importance weights collapse (lesson 14).
• "What is a good acceptance rate?" — About 0.44 for 1-D random-walk Metropolis, about 0.234 in high dimension, about 0.65 to 0.8 for HMC. A rate near 1 means tiny steps.
• "Your chains have $\hat R = 1.2$. What do you do?" — Look at the trace plots, run longer, reparametrise (for example a non-centred parametrisation of a hierarchical model), use a better sampler, and check again. Do not drop chains.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. stationary_distribution by a linear solve and by an eigenvector of $P^\top$ (not by $P^t$), tested on the weather chain, the periodic flip chain and a birth–death chain.
2. metropolis_hastings with the log acceptance test; the acceptance rate against $\frac2\pi\arctan(2/s)$ and on a banana-shaped target.
3. gibbs_bivariate_gaussian and the lag-1 autocorrelation $\rho^2$ of its chain, for several $\rho$.
4. leapfrog, tested for reversibility and the $O(\varepsilon^2)$ energy error; then a small HMC sampler on a correlated Gaussian.
5. ess_autocorr (Geyer) and split_rhat, tested on AR(1) chains ($\text{ESS}/n = (1 - \varphi)/(1 + \varphi)$) and on chains that are stuck or drift.

11. Watch and read #

Lecture 31: Markov Chains | Statistics 110Harvard University · Joe Blitzstein, Statistics 110 · 47 min

The Markov property, transition matrices, $\mu_{t+1} = \mu_t P$ and the stationary distribution as an eigenvector.

Lecture 32: Markov Chains Continued | Statistics 110Harvard University · Joe Blitzstein, Statistics 110 · 48 min

Irreducible and periodic chains, reversibility (detailed balance) and random walks on graphs.

16. Markov Chains IMIT OpenCourseWare · John Tsitsiklis, MIT 6.041 · 52 min

A second, slower introduction with many small examples of states, transitions and n-step probabilities.

Markov Chain Monte Carlo (MCMC) : Data Science Conceptsritvikmath · 12 min

Why we build a Markov chain whose stationary distribution is the target, and the role of burn-in.

Metropolis - Hastings : Data Science Conceptsritvikmath · 18 min

The acceptance probability, and why it makes the target the stationary distribution of the chain.

An introduction to Gibbs samplingBen Lambert · 19 min

How Gibbs sampling updates one coordinate at a time from its conditional distribution.

The intuition behind the Hamiltonian Monte Carlo algorithmBen Lambert · 32 min

The physical picture behind HMC, the role of the momentum, and why it explores faster than a random walk.

• Joseph Blitzstein and Jessica Hwang, Introduction to Probability (free online, 2nd ed.): ch. 11 (Markov chains: classification of states, stationary distributions, reversibility) and ch. 12 (Markov chain Monte Carlo: Metropolis–Hastings and Gibbs).
• Art Owen, Monte Carlo theory, methods and examples (free chapters): ch. 11 (Markov chains, detailed balance, Metropolis–Hastings, burn-in and error estimates) and ch. 12 (Gibbs sampling and correlated components).
• Michael Betancourt, A Conceptual Introduction to Hamiltonian Monte Carlo (2017): §3 (phase space and Hamilton's equations) and §5 (symplectic integrators, the Metropolis correction, divergences).
• Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter and Paul-Christian Bürkner, Rank-normalization, folding, and localization: an improved $\hat R$ (2021): §2 (recommendations), §3 (split-$\hat R$ and ESS with Geyer's truncation).
• Radford Neal, MCMC using Hamiltonian dynamics (2011): the classic tutorial on HMC, with the leapfrog proofs and tuning advice.