Markov Chain Monte Carlo (MCMC): Metropolis-Hastings and Gibbs Sampling

Introduction

In Bayesian inference, sampling from the posterior distribution \(p(\theta | D)\) is a central challenge. In most cases, this distribution is analytically intractable and cannot be directly sampled.

Markov Chain Monte Carlo (MCMC) constructs a Markov chain whose stationary distribution is the target distribution, generating samples from it. Running the chain long enough yields approximate samples from the posterior.

Markov Chain Basics

A Markov chain is a stochastic process where the next state depends only on the current state:

The goal of MCMC is to design a Markov chain with target distribution \(\pi(\theta)\) as its stationary distribution, guaranteed by the detailed balance condition:

Why Detailed Balance Guarantees Stationarity

Saying a Markov chain has \(\pi\) as its stationary distribution means the following relation holds:

The left-hand side is the marginal distribution of the state \(\theta'\) after one transition, assuming the current state is distributed according to \(\pi\) . If this matches the original \(\pi(\theta')\) , the distribution’s shape is unchanged by the transition.

Integrating both sides of the detailed balance condition (2) over \(\theta\) gives

which is exactly equation (2’). The last equality uses the fact that \(T(\theta | \theta')\) is a transition kernel starting from \(\theta'\) , so integrating over the destination \(\theta\) must equal 1 (normalization of a transition probability). Hence any kernel satisfying detailed balance automatically satisfies stationarity.

The converse does not hold: a stationary distribution need not satisfy detailed balance (reversibility). However, reversibility is a sufficient condition for stationarity and is tractable to design transition kernels around, which is why the major MCMC methods (Metropolis-Hastings, Gibbs sampling, and HMC discussed in a later article) all guarantee stationarity by construction through reversibility.

Metropolis-Hastings Algorithm

Algorithm

1. Set initial value \(\theta_0\)
2. Generate candidate \(\theta'\) from proposal distribution \(q(\theta' | \theta_t)\)
3. Compute acceptance probability:

4. Accept \(\theta_{t+1} = \theta'\) with probability \(\alpha\) ; otherwise \(\theta_{t+1} = \theta_t\)
5. Repeat steps 2-4

For symmetric proposals (\(q(\theta' | \theta) = q(\theta | \theta')\) ), this simplifies to the Metropolis algorithm:

Crucially, the normalization constant of \(\pi(\theta)\) is not needed. In Bayesian inference, only the unnormalized density \(p(\theta | D) \propto p(D | \theta) p(\theta)\) is required.

Proof That the Acceptance Probability Satisfies Detailed Balance

Rather than accepting equation (3) at face value, let’s prove why the acceptance probability must take this exact form. The MH transition kernel has a component where a proposal \(\theta' \neq \theta\) is accepted, and a component where it is rejected and the chain stays put. For \(\theta' \neq \theta\) ,

So it suffices to show detailed balance \(\pi(\theta) T(\theta' | \theta) = \pi(\theta') T(\theta | \theta')\) holds for \(\theta' \neq \theta\) , i.e.,

Define the ratio \(r = \dfrac{\pi(\theta') q(\theta | \theta')}{\pi(\theta) q(\theta' | \theta)}\) , and set \(\alpha(\theta' | \theta) = \min(1, r)\) , \(\alpha(\theta | \theta') = \min(1, 1/r)\) as in equation (3).

Case \(r \le 1\) : \(\alpha(\theta' | \theta) = r\) and \(\alpha(\theta | \theta') = 1\) , so

Case \(r > 1\) : \(\alpha(\theta' | \theta) = 1\) , \(\alpha(\theta | \theta') = 1/r\) , and the symmetric calculation shows both sides equal \(\pi(\theta) q(\theta' | \theta)\) . When \(\theta' = \theta\) , both sides trivially contain the identical “stay put” term.

Detailed balance therefore holds in every case. This shows equation (3) is not an arbitrarily chosen formula, but the design that satisfies detailed balance while minimizing rejections (keeping \(\alpha\) as close to 1 as the constraint allows).

Python Implementation

Sampling from a mixture of Gaussians:

import numpy as np
import matplotlib.pyplot as plt

def metropolis_hastings(log_target, initial, n_samples, proposal_std=1.0):
    """Metropolis-Hastings algorithm"""
    samples = [initial]
    current = initial
    accepted = 0

    for _ in range(n_samples):
        proposal = current + np.random.normal(0, proposal_std)
        log_alpha = log_target(proposal) - log_target(current)

        if np.log(np.random.random()) < log_alpha:
            current = proposal
            accepted += 1

        samples.append(current)

    acceptance_rate = accepted / n_samples
    return np.array(samples), acceptance_rate

# --- Target: Gaussian mixture ---
def log_target(x):
    """Log target density: 0.3*N(-2,1) + 0.7*N(3,0.5)"""
    from scipy.special import logsumexp
    log_p1 = np.log(0.3) - 0.5 * (x + 2)**2
    log_p2 = np.log(0.7) - (x - 3)**2
    return logsumexp([log_p1, log_p2])

# --- Run ---
np.random.seed(42)
samples, acc_rate = metropolis_hastings(log_target, initial=0.0,
                                        n_samples=50000, proposal_std=1.5)
print(f"Acceptance rate: {acc_rate:.3f}")

burn_in = 5000
samples = samples[burn_in:]

# --- Visualization ---
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].plot(samples[:2000], alpha=0.5, linewidth=0.5)
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('Sample value')
axes[0].set_title('Trace Plot')
axes[0].grid(True, alpha=0.3)

x = np.linspace(-6, 6, 200)
true_pdf = 0.3 * np.exp(-0.5 * (x + 2)**2) / np.sqrt(2*np.pi) + \
           0.7 * np.exp(-0.5 * ((x - 3)/0.5)**2) / (0.5 * np.sqrt(2*np.pi))
axes[1].hist(samples, bins=100, density=True, alpha=0.5, label='MCMC samples')
axes[1].plot(x, true_pdf, 'r-', linewidth=2, label='True density')
axes[1].set_xlabel('x')
axes[1].set_ylabel('Density')
axes[1].set_title('Histogram vs True Distribution')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Gibbs Sampling

For multidimensional cases, sample each variable sequentially from its conditional distribution. For two variables \((\theta_1, \theta_2)\) :

1. \(\theta_1^{(t+1)} \sim p(\theta_1 | \theta_2^{(t)}, D)\)
2. \(\theta_2^{(t+1)} \sim p(\theta_2 | \theta_1^{(t+1)}, D)\)

Particularly effective when conditional distributions are known (e.g., conjugate priors). It is a special case of Metropolis-Hastings with acceptance rate always 1.

Derivation as a Special Case of Metropolis-Hastings

The fact that the acceptance rate is always 1 is not an assumption to accept at face value — it follows directly from the MH framework. For a joint distribution \(\pi(\theta_1, \theta_2)\) , take the proposal for updating \(\theta_1\) to be the full conditional distribution itself:

(with \(\theta_2\) held fixed, \(\theta_1\) is sampled directly from its conditional distribution. The reverse proposal \(q(\theta_1 | \theta_1', \theta_2)\) follows the same conditional \(\pi(\theta_1 | \theta_2)\) .)

The joint distribution factors as \(\pi(\theta_1, \theta_2) = \pi(\theta_1 | \theta_2) \pi(\theta_2)\) . Substituting into the MH acceptance probability (3):

The factors \(\pi(\theta_1' | \theta_2)\) , \(\pi(\theta_1 | \theta_2)\) , and \(\pi(\theta_2)\) in the numerator and denominator all cancel, so the ratio is identically 1. This proves that Gibbs sampling proposals are always accepted. Sampling directly from the conditional distribution is itself already the choice of proposal that maximizes the MH acceptance probability.

Ergodicity and the Theory of Convergence

Detailed balance guarantees the existence of a stationary distribution, but on its own it does not guarantee that the chain will actually converge to it. Convergence requires two additional properties.

Irreducibility

From any state \(\theta\) , the chain must be able to reach any state \(\theta'\) with \(\pi(\theta') > 0\) in a finite number of steps with positive probability. Formally, there must exist some \(n\) such that the \(n\) -step transition kernel \(T^{(n)}\) satisfies

If irreducibility fails, the chain can become permanently confined to part of the state space, never reaching some of the target distribution’s mass.

Aperiodicity

The greatest common divisor of the set of possible return times to a given state must be 1. A chain with period \(d > 1\) cycles among \(d\) distinct distributions at time \(t\) and never converges to a single limiting distribution (this can happen, for instance, if the set of states visited at even versus odd steps is designed to be mutually exclusive). As long as a random-walk proposal (additive noise from a continuous distribution) is used, aperiodicity is typically satisfied automatically.

The Ergodic Theorem

If a chain is irreducible and aperiodic with stationary distribution \(\pi\) , then starting from any initial distribution, the distribution at time \(n\) converges to \(\pi\) in total variation distance:

Moreover, as an extension of the law of large numbers, for any integrable function \(f\) the time average converges almost surely to the true expectation (the ergodic average):

This is the theoretical justification at the heart of MCMC. But the theorem only guarantees eventual convergence — it says nothing about how fast (the mixing time). The two concrete examples below show why this distinction matters in practice.

Example 1: A Case Where Irreducibility Is Strictly Violated

Consider a target distribution that is uniform over two disjoint intervals, \(\pi(\theta) = 0.5 \cdot \mathrm{Unif}[0, 1] + 0.5 \cdot \mathrm{Unif}[3, 4]\) , with a bounded uniform random-walk proposal \(q(\theta' | \theta) = \mathrm{Unif}(\theta - 0.4, \theta + 0.4)\) (a single step moves by at most 0.4). From any point in \([0, 1]\) , the proposal’s reach is at most \((-0.4, 1.4)\) , which never overlaps \([3, 4]\) . This is not an approximation issue — \(T^{(n)}(\theta' | \theta) = 0\) holds exactly for every \(n\) (with \(\theta \in [0,1]\) , \(\theta' \in [3,4]\) ), a complete violation of irreducibility. The chain is permanently confined to whichever interval it starts in, entirely missing the other 50% of the target’s mass. No simulation is needed here — it follows immediately from the geometry of the proposal’s reach.

Example 2: “Practical” Trapping in Real Runs

The example above is a case where irreducibility fails strictly. What comes up far more often in practice is a chain that is theoretically irreducible (a Gaussian proposal has positive density everywhere on \(\mathbb{R}\) , so every state is reachable with positive probability), yet is effectively trapped within any realistic computation time because the step size is too small relative to the width of the valley in the target distribution.

To verify this concretely, we used the same mixture-of-Gaussians target from the Python implementation above, \(0.3 \cdot \mathcal{N}(-2, 1) + 0.7 \cdot \mathcal{N}(3, 0.5)\) , starting at \(\theta_0 = -2\) (the minor mode) and running for 2000 iterations at two step sizes (proposal standard deviations), 0.1 and 1.5, with the random seed and all other settings held fixed.

At step size 0.1, the acceptance rate is a seemingly “efficient” 97.0%, yet the chain never once reaches the dominant mode (near \(x=3\) ) over 2000 iterations, remaining completely trapped near the minor mode (sample mean \(-1.97\) ). At step size 1.5, the chain crosses between modes 48 times and explores both.

Tracking the convergence of the fraction with \(x > 0.5\) as we extend the run length (step size 1.5, same start at \(-2\) ; the true value of this fraction is about 0.625, as explained below):

Even with a “good” step size, a multimodal target can take on the order of hundreds of thousands to millions of iterations to converge to the true mixing proportions. This is a concrete illustration of the gap between what the ergodic theorem guarantees (“eventual” convergence) and what is achievable within a practically acceptable amount of computation time.

Implementation caveat (a normalization-constant pitfall): The log_target function in this article’s Python code implements only the exponent of each component, \(\log p_1 = \log(0.3) - 0.5(x+2)^2\) and \(\log p_2 = \log(0.7) - (x-3)^2\) , omitting each Gaussian’s proper normalization constant (\(1/\sqrt{2\pi}\) for \(\mathcal{N}(-2,1)\) , \(1/\sqrt{\pi}\) for \(\mathcal{N}(3, 0.5)\) ). MH does not require knowing the normalization constant shared by the entire target density, but the relative normalization between mixture components directly affects the resulting mixing proportions. Numerically integrating this density shows that, despite the “0.3” and “0.7” labels in the code, the true mixing proportions are approximately \(0.377 : 0.623\) (giving \(P(x > 0.5) \approx 0.625\) ). This is exactly the value the experimental fractions above are converging toward. When implementing a mixture density, always include each component’s normalization constant — a classic pitfall otherwise.

Convergence Diagnostics

The ergodic theorem guarantees eventual convergence, but a separate, quantitative way to judge whether a finite chain in hand has “converged enough” is still needed. Below are two standard diagnostics, along with numbers computed by actually running the Python code.

The Gelman-Rubin Statistic \(\hat{R}\)

Run \(m\) independent chains from different starting points, each of length \(n\) , and compare between-chain and within-chain variance. Writing the mean of chain \(j\) as \(\bar{\theta}_j\) and the grand mean as \(\bar{\theta}\) , the average within-chain variance is

and the between-chain variance (scaled by \(n\) to match the scale of the within-chain variance) is

Using these to form the pooled variance estimate \(\widehat{\mathrm{Var}}(\theta) = \frac{n-1}{n} W + \frac{1}{n} B\) , define

If every chain has converged to the stationary distribution, each chain mean \(\bar{\theta}_j\) approaches the same mean of \(\pi\) , so \(B\) shrinks and \(\hat{R} \to 1\) . Conversely, if the chains still carry the imprint of their different starting points, the \(\bar{\theta}_j\) disagree, \(B\) grows, and \(\hat{R} \gg 1\) . A common rule of thumb treats \(\hat{R} < 1.1\) (or \(1.01\) under a stricter standard) as evidence of convergence.

Experimental verification: For the mixture-of-Gaussians target, we ran 4 chains (step size 1.5) from widely separated starting points \(\theta_0 \in \{-10, -5, 5, 10\}\) and computed \(\hat{R}\) as a function of chain length \(n\) (with no burn-in removed).

At \(n=50\) –\(100\) , the differing starting points still dominate, giving a very large \(\hat{R}\) of 1.85–1.89. By \(n=300\) it drops sharply to 1.10, and from \(n=1000\) onward it stays below 1.02. This is consistent with the burn-in results below, which show that the influence of extreme starting points dissipates within tens to hundreds of iterations.

Fixing the chain length at \(n=2000\) and instead varying the amount of burn-in discarded:

Once \(n\) reaches 2000, \(\hat{R}\) is already below 1.01 even with no burn-in removed at all — in this setting, a sufficiently long chain already dilutes the effect of initial-value dependence on its own. (This is specific to an example where each chain reaches a reasonable region within a few dozen iterations; removing burn-in remains important in general for shorter chains.)

Effective Sample Size (ESS)

Because MCMC samples are autocorrelated, \(n\) samples do not carry as much information as \(n\) independent samples. The effective sample size (ESS) is the “equivalent number of independent samples” that would give the same estimation precision as the actual autocorrelated sample. Writing the lag-\(k\) autocorrelation as \(\rho_k\) , define the integrated autocorrelation time

so that \(\mathrm{ESS} = n / \tau_{\mathrm{int}}\) . In practice the infinite sum must be truncated; this article uses Geyer’s initial positive sequence rule (summing pairs \(\rho_{2k} + \rho_{2k+1}\) as long as the pair remains positive, then stopping). The next section uses this ESS as the metric for choosing step size.

Practical Considerations

Burn-in

Discard the first several thousand samples (burn-in period) to remove the influence of the initial value. To verify this need directly, we ran chains on the mixture-of-Gaussians target from 5 widely separated initial values, \(\theta_0 \in \{-20, -10, 0, 10, 20\}\) (the target’s mass lies mostly in roughly \([-4, 5]\) ), and measured the number of iterations before each chain first entered that plausible region (step size 1.5).

Every initial value reaches the plausible region within 40 iterations, confirming that the burn-in period of 5000 samples used in this article’s Python implementation is a very conservative choice for a target this low-dimensional and simple. That said, for high-dimensional, strongly multimodal, or highly correlated posteriors, reaching a plausible region may take only a few dozen iterations while the chain’s statistics (variance, higher moments) can still take far longer to converge — the \(\hat{R}\) and ESS diagnostics from the previous section quantify “degree of statistical convergence,” not merely “has it reached a plausible region,” and the two should not be conflated.

Acceptance Rate Tuning

Too large a proposal standard deviation lowers the acceptance rate; too small slows exploration. To quantify this trade-off, we used a standard normal \(\mathcal{N}(0,1)\) target, swept the step size (proposal standard deviation) from 0.1 to 10, generated 50,000 samples, discarded a burn-in of 5,000, and computed the acceptance rate and ESS on the remaining 45,000 samples.

ESS is an inverted-U function of step size, peaking at 10,508.6 (23.35% of the sample count) at step size 2.38 (acceptance rate 44.8%). At too small a step size (0.1), nearly every proposal is accepted, but each step moves so little that strong autocorrelation remains (\(\tau_{\mathrm{int}} \approx 469\) ), dropping ESS to about 0.2% of the sample count. At too large a step size (10.0), most proposals land in low-probability regions and are rejected, dropping the acceptance rate to 12.4% and ESS falls again.

The acceptance rate of 44.8% at this optimum agrees closely with the theoretical optimum of about 44% derived by Roberts, Gelman & Gilks (1997) for a one-dimensional i.i.d. target (in high dimensions with strong correlation, this optimum is known to approach about 23.4% as dimensionality grows). The practical rule of thumb of tuning the proposal’s variance “neither too large nor too small” is reproduced directly by this numerical experiment.

Autocorrelation and Thinning

Consecutive samples are correlated. When independent samples are needed, apply thinning (keeping every \(k\) -th sample). However, the ESS defined above already quantifies the information needed to correctly estimate statistics such as variance without thinning, so in practice it is more efficient to size the run based on ESS rather than to thin. Thinning is useful for reducing storage, but since the information in discarded samples is lost entirely, keeping every sample is (for the same computational cost) more efficient for variance estimation.

Related Articles

• Hamiltonian Monte Carlo (HMC) in Python - A follow-up solving the autocorrelation problem noted here using gradient information and Hamiltonian dynamics, measuring roughly a 45x efficiency gain on a correlated Gaussian.
• Simulated Annealing: Theory and Python Implementation - Applies the same Metropolis criterion to optimization.
• Bayesian Optimization: Fundamentals and Python Implementation - Using MCMC-estimated posteriors in Bayesian optimization.
• Particle Filter Python Implementation: Comparing Resampling Methods - Sequential Monte Carlo (SMC), closely related to MCMC.
• Bayesian Linear Regression: From Least Squares to Bayesian Estimation - Analytical Bayesian inference without MCMC.
• Cross-Entropy Method: A Practical Monte Carlo Optimization Technique - Monte Carlo methods applied to optimization.
• Gaussian Process Regression: Theory and Python Implementation - A flagship Bayesian model; contrasts MCMC’s general sampling approach with the GP’s closed-form posterior.

References

• Metropolis, N., et al. (1953). “Equation of State Calculations by Fast Computing Machines”. The Journal of Chemical Physics, 21(6), 1087-1092.
• Hastings, W. K. (1970). “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”. Biometrika, 57(1), 97-109.
• Gelman, A., et al. (2013). Bayesian Data Analysis (3rd ed.). Chapman and Hall/CRC. Chapters 11-12.
• Gelman, A., & Rubin, D. B. (1992). “Inference from Iterative Simulation Using Multiple Sequences”. Statistical Science, 7(4), 457-472.
• Roberts, G. O., Gelman, A., & Gilks, W. R. (1997). “Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms”. The Annals of Applied Probability, 7(1), 110-120.
• Geyer, C. J. (1992). “Practical Markov Chain Monte Carlo”. Statistical Science, 7(4), 473-483.