Markov chain Monte Carlo: derivation, tuning, and diagnostics

NoteWhat you should be able to do after this chapter
  • explain why MCMC is needed for posterior integration;
  • derive the Metropolis-Hastings acceptance probability;
  • describe Gibbs sampling as exact conditional updating;
  • tune a random-walk proposal intelligently rather than heuristically;
  • interpret autocorrelation, ESS, trace plots, and \(\hat{R}\).

1 The computational problem

In Bayesian inference, the target quantity is often an expectation under the posterior:

\[ \mathbb{E}[h(\theta) \mid y] = \int h(\theta)\, p(\theta \mid y)\, d\theta. \]

Outside conjugate or low-dimensional settings, this integral is unavailable in closed form. MCMC solves the problem by constructing a Markov chain whose stationary distribution is the posterior (Robert and Casella 2004; Brooks et al. 2011).

If \(\theta^{(1)}, \dots, \theta^{(S)}\) are approximately distributed from \(p(\theta \mid y)\), then

\[ \mathbb{E}[h(\theta) \mid y] \approx \frac{1}{S}\sum_{s=1}^S h(\theta^{(s)}). \]

The subtlety is that the draws are not independent. Dependence affects Monte Carlo error, which is why diagnostics are central rather than optional.

2 Markov chains and invariance

A Markov chain with transition kernel \(K(\theta, d\theta')\) is a stochastic process in which

\[ \Pr(\theta_{t+1} \in A \mid \theta_t, \theta_{t-1}, \dots) = \Pr(\theta_{t+1} \in A \mid \theta_t). \]

We want a chain whose invariant distribution is the posterior \(\pi(\theta)=p(\theta \mid y)\), meaning

\[ \int \pi(\theta) K(\theta, d\theta')\, d\theta = \pi(\theta'). \]

One sufficient condition is detailed balance.

Proposition 1 If

\[ \pi(\theta) K(\theta, d\theta') = \pi(\theta') K(\theta', d\theta) \]

for all states, then \(\pi\) is invariant for the Markov chain.

2.1 Proof sketch

Integrating both sides over \(\theta\) gives

\[ \int \pi(\theta) K(\theta, d\theta')\, d\theta = \int \pi(\theta') K(\theta', d\theta)\, d\theta = \pi(\theta'), \]

because the kernel integrates to one in its second argument.

Detailed balance is stronger than necessary, but it provides an elegant route to MCMC design.

3 Metropolis and Metropolis-Hastings

Suppose the posterior is known up to a normalizing constant:

\[ \pi(\theta) \propto \tilde{\pi}(\theta). \]

Given the current state \(\theta\), propose \(\theta' \sim q(\theta' \mid \theta)\) and accept it with probability

\[ \alpha(\theta, \theta') = \min\left\{1, \frac{\tilde{\pi}(\theta') q(\theta \mid \theta')} {\tilde{\pi}(\theta) q(\theta' \mid \theta)} \right\}. \]

If the proposal is rejected, stay at \(\theta\).

3.1 Derivation of the acceptance ratio

The transition kernel has two components:

  • a move from \(\theta\) to \(\theta'\) with density \(q(\theta' \mid \theta)\alpha(\theta,\theta')\);
  • a point mass at \(\theta\) representing rejection.

For \(\theta \neq \theta'\), detailed balance requires

\[ \pi(\theta) q(\theta' \mid \theta)\alpha(\theta,\theta') = \pi(\theta') q(\theta \mid \theta')\alpha(\theta',\theta). \]

Let

\[ r(\theta,\theta') = \frac{\pi(\theta') q(\theta \mid \theta')} {\pi(\theta) q(\theta' \mid \theta)}. \]

Choose

\[ \alpha(\theta,\theta')=\min\{1,r(\theta,\theta')\}. \]

Then:

  • if \(r \le 1\), the left-hand side becomes \(\pi(\theta) q(\theta' \mid \theta)r = \pi(\theta') q(\theta \mid \theta')\) and the reverse acceptance is \(1\);
  • if \(r > 1\), the left-hand side is \(\pi(\theta) q(\theta' \mid \theta)\) and the right-hand side is \(\pi(\theta') q(\theta \mid \theta') (1/r)\).

Thus detailed balance holds in both cases (Hastings 1970).

3.2 Symmetric proposals

If

\[ q(\theta' \mid \theta)=q(\theta \mid \theta'), \]

then the proposal terms cancel and

\[ \alpha(\theta,\theta') = \min\left\{1, \frac{\tilde{\pi}(\theta')}{\tilde{\pi}(\theta)}\right\}. \]

This is the usual random-walk Metropolis rule.

4 Gibbs sampling

Gibbs sampling updates components from their exact full conditional distributions. If \(\theta=(\theta_1,\dots,\theta_p)\), a one-at-a-time Gibbs sweep draws

\[ \theta_1^{(t+1)} \sim p(\theta_1 \mid \theta_2^{(t)},\dots,\theta_p^{(t)}, y), \]

then

\[ \theta_2^{(t+1)} \sim p(\theta_2 \mid \theta_1^{(t+1)},\theta_3^{(t)},\dots,\theta_p^{(t)}, y), \]

and so on.

Because each step draws from the correct conditional, Gibbs sampling inherits the posterior as its invariant distribution (Geman and Geman 1984).

Gibbs is powerful when full conditionals are standard. It is limited when they are not.

5 Example: random-walk Metropolis for a Normal posterior

Consider a simple model with known variance:

\[ y_i \mid \theta \overset{\text{i.i.d.}}{\sim} \mathcal{N}(\theta, \sigma^2), \qquad \theta \sim \mathcal{N}(\mu_0, \tau_0^2). \]

The posterior is Gaussian:

\[ \theta \mid y \sim \mathcal{N}(\mu_n, \tau_n^2), \]

where

\[ \tau_n^2 = \left(\frac{n}{\sigma^2}+\frac{1}{\tau_0^2}\right)^{-1}, \qquad \mu_n = \tau_n^2 \left(\frac{n\bar{y}}{\sigma^2}+\frac{\mu_0}{\tau_0^2}\right). \]

This model is useful pedagogically because we know the exact answer and can therefore diagnose MCMC performance honestly.

5.1 Unnormalized log posterior

Ignoring additive constants,

\[ \log \tilde{\pi}(\theta) = -\frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-\theta)^2 -\frac{1}{2\tau_0^2}(\theta-\mu_0)^2. \]

5.2 R implementation

The following code implements a random-walk Metropolis sampler. A reusable script is also available at mh_sampler.R.

log_post_theta <- function(theta, y, sigma, mu0, tau0) {
  loglik <- sum(dnorm(y, mean = theta, sd = sigma, log = TRUE))
  logprior <- dnorm(theta, mean = mu0, sd = tau0, log = TRUE)
  loglik + logprior
}

rw_metropolis <- function(init, n_iter, proposal_sd, y, sigma, mu0, tau0) {
  draws <- numeric(n_iter)
  draws[1] <- init
  accept <- 0L

  for (t in 2:n_iter) {
    proposal <- rnorm(1, mean = draws[t - 1], sd = proposal_sd)
    log_alpha <- log_post_theta(proposal, y, sigma, mu0, tau0) -
      log_post_theta(draws[t - 1], y, sigma, mu0, tau0)

    if (log(runif(1)) < log_alpha) {
      draws[t] <- proposal
      accept <- accept + 1L
    } else {
      draws[t] <- draws[t - 1]
    }
  }

  list(draws = draws, accept_rate = accept / (n_iter - 1))
}

set.seed(42)
y <- rnorm(40, mean = 1.5, sd = 2)
fit <- rw_metropolis(
  init = 0,
  n_iter = 10000,
  proposal_sd = 0.35,
  y = y,
  sigma = 2,
  mu0 = 0,
  tau0 = 5
)

fit$accept_rate

6 Tuning proposal scale

Random-walk Metropolis reveals a basic MCMC trade-off:

  • if the proposal scale is too small, acceptance is high but exploration is slow;
  • if the proposal scale is too large, proposals often land in low-density regions and are rejected.

Good MCMC is not about maximizing acceptance rate. It is about maximizing effective exploration per unit computation.

In one dimension, an acceptance rate around \(0.4\) to \(0.5\) is often reasonable for random walks; in higher dimensions, lower optimal rates can occur. These are rules of thumb, not universal laws.

7 Autocorrelation and effective sample size

Because MCMC draws are dependent, the Monte Carlo error of the sample mean depends on autocorrelation. If \(\rho_k\) denotes lag-\(k\) autocorrelation, then asymptotically

\[ \operatorname{Var}\!\left(\frac{1}{S}\sum_{s=1}^S h(\theta^{(s)})\right) \approx \frac{\sigma_h^2}{S}\left(1 + 2\sum_{k=1}^{\infty}\rho_k\right). \]

This motivates the effective sample size

\[ \operatorname{ESS} = \frac{S}{1 + 2\sum_{k=1}^{\infty}\rho_k}. \]

If autocorrelation is high, ESS can be much smaller than the nominal number of iterations.

8 Diagnostics that actually matter

8.1 Trace plots

Useful for detecting:

  • nonstationarity;
  • sticking;
  • multimodality or poor mixing;
  • obvious chain-to-chain disagreement.

Not sufficient for certifying convergence.

8.2 Acceptance rate

Useful for diagnosing an obviously bad proposal scale. Not sufficient for certifying efficiency or correctness.

8.3 \(\hat{R}\)

Multiple chains should mix well both within and across chains. Modern rank-normalized \(\hat{R}\) near 1 is a necessary but not sufficient condition for trustworthy inference (Vehtari et al. 2021).

8.4 ESS

ESS is a measure of Monte Carlo information, not a measure of model adequacy.

8.5 Domain-specific diagnostics

For HMC, divergences and E-BFMI matter. For importance sampling, tail weights matter. Diagnostics must match the algorithm.

9 Common failure modes

9.1 Poor scaling

If parameters are on very different scales, a simple random-walk proposal can mix terribly.

9.2 Multi-modality

A chain can appear stable while exploring only one mode.

9.3 Strong posterior correlation

Random-walk proposals become inefficient on narrow, tilted ridges. This is a preview of why HMC is valuable.

9.4 False reassurance from long runs

A very long chain with high autocorrelation can still carry little effective information.

10 Gibbs sampling example sketch

Suppose

\[ y_i \mid \mu,\tau \sim \mathcal{N}(\mu,\tau^{-1}), \qquad \mu \mid \tau \sim \mathcal{N}(\mu_0,(\kappa_0\tau)^{-1}), \qquad \tau \sim \operatorname{Gamma}(a_0,b_0). \]

Then the full conditionals are:

  • \(\mu \mid \tau, y\) is normal;
  • \(\tau \mid \mu, y\) is gamma.

This produces a two-step Gibbs sampler. The point is not the algebra alone, but the structural insight: conjugacy turns conditionals into direct samplers, which is why Gibbs flourished historically.

11 Exercises

  1. Derive the Metropolis-Hastings acceptance ratio from detailed balance.
  2. Show that Gibbs sampling leaves the posterior invariant.
  3. Simulate a random-walk Metropolis chain for the normal posterior above with three proposal scales. Compare acceptance rate, ESS, and empirical autocorrelation.
  4. Construct an independence sampler for a unimodal target and compare it with a random walk.
  5. Explain why thinning is usually a poor substitute for running an efficient sampler.

12 Summary

MCMC works by designing a Markov chain whose invariant distribution is the posterior. The most important ideas in this chapter are:

  • invariance through detailed balance;
  • Metropolis-Hastings as accept-reject correction of a proposal mechanism;
  • Gibbs sampling as exact conditional replacement;
  • Monte Carlo error under dependence;
  • diagnostics as essential components of inference rather than optional afterthoughts.

The next computational step is to move beyond random walks and exploit posterior geometry, which is the subject of advanced MCMC and Stan workflow.

13 Further reading

The original papers are short and still readable. Metropolis et al. (1953) introduced the symmetric-proposal algorithm for a physics problem; Hastings (1970) generalised it to asymmetric proposals and gave the acceptance ratio in the form used today; Geman and Geman (1984) introduced Gibbs sampling for image restoration, which is where the name comes from.

Robert and Casella (2004) is the standard textbook treatment and covers the convergence theory this chapter states without proof — irreducibility, aperiodicity, and the ergodic theorem for Markov chains on general state spaces. Brooks et al. (2011) is the handbook, with survey chapters on most of the variations you are likely to meet.

Andrieu et al. (2003) is a good introduction for readers coming from machine learning. On diagnostics, Vehtari et al. (2021) is the current reference for rank-normalised \(\widehat{R}\) and effective sample size, and supersedes the older diagnostics still found in some software. The geometric methods that overcome the random-walk limitations described here are covered in advanced MCMC, and the mixture-specific case where the dimension itself varies in the telescoping sampler.

Back to top

References

Andrieu, Christophe, Nando de Freitas, Arnaud Doucet, and Michael I. Jordan. 2003. “An Introduction to MCMC for Machine Learning.” Machine Learning 50 (1–2): 5–43.
Brooks, Steve, Andrew Gelman, Galin Jones, and Xiao-Li Meng. 2011. Handbook of Markov Chain Monte Carlo. CRC Press.
Geman, Stuart, and Donald Geman. 1984. “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images.” IEEE Transactions on Pattern Analysis and Machine Intelligence 6 (6): 721–41.
Hastings, W. K. 1970. “Monte Carlo Sampling Methods Using Markov Chains and Their Applications.” Biometrika 57 (1): 97–109.
Metropolis, Nicholas, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. 1953. “Equation of State Calculations by Fast Computing Machines.” Journal of Chemical Physics 21 (6): 1087–92.
Robert, Christian P., and George Casella. 2004. Monte Carlo Statistical Methods. 2nd ed. Springer.
Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2021. “Rank-Normalization, Folding, and Localization: An Improved \(\hat{R}\) for Assessing Convergence of MCMC.” Bayesian Analysis 16 (2): 667–718.