MCMC Diagnostics and Convergence

Markov Chain Monte Carlo (MCMC) methods generate samples from posterior distributions that cannot be evaluated analytically. Proper diagnostics are essential to ensure the samples are valid representations of the target distribution.

MCMC Fundamentals

Metropolis-Hastings Algorithm

Gibbs Sampling

Convergence Diagnostics

R-hat (Gelman-Rubin Statistic)

More precisely:

Effective Sample Size (ESS)

ESS measures the number of independent samples equivalent to the correlated MCMC output. High autocorrelation reduces ESS. A minimum of ESS is recommended for reliable posterior summaries.

Autocorrelation

Thinning and Burn-in

Geweke Diagnostic

Brooks-Gelman-Rubin Multivariate Diagnostic

The multivariate (Vehtari et al., 2021) monitors all parameters simultaneously and is more sensitive to divergence in high-dimensional posteriors. is the recommended threshold.

Divergent Transitions (HMC/NUTS)

Python Implementation

import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Simple model for demonstration ---
n_obs = 200
X = np.random.randn(n_obs)
true_beta = np.array([1.5, -2.0])
y = true_beta[0] + true_beta[1] * X + np.random.randn(n_obs) * 0.5

with pm.Model() as demo_model:
    beta = pm.Normal('beta', mu=0, sigma=10, shape=2)
    sigma = pm.HalfNormal('sigma', sigma=2)
    mu = beta[0] + beta[1] * X
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
    
    # Run MCMC with diagnostics
    trace = pm.sample(2000, tune=1000, chains=4, return_inferencedata=True,
                      random_seed=42)

# --- R-hat ---
print("=== R-hat ===")
print(az.rhat(trace, var_names=['beta', 'sigma']))

# --- Effective Sample Size ---
print("\n=== Effective Sample Size ===")
print(az.ess(trace, var_names=['beta', 'sigma']))

# --- Autocorrelation ---
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for i, ax in enumerate(axes.flatten()):
    if i < 2:
        az.plot_autocorr(trace, var_names=['beta'], combined=True, ax=ax)
    elif i == 2:
        az.plot_autocorr(trace, var_names=['sigma'], combined=True, ax=ax)
    else:
        az.plot_posterior(trace, var_names=['beta'], ax=ax)
plt.tight_layout()
plt.savefig('mcmc_autocorr.png', dpi=150)
plt.show()

# --- Trace plots ---
az.plot_trace(trace, var_names=['beta', 'sigma'], compact=True)
plt.tight_layout()
plt.savefig('mcmc_trace.png', dpi=150)
plt.show()

# --- Rank plots (modern diagnostic) ---
az.plot_rank(trace, var_names=['beta'])
plt.tight_layout()
plt.savefig('mcmc_rank.png', dpi=150)
plt.show()

# --- Geweke diagnostic ---
# Manual implementation for demonstration
from scipy import stats as sp_stats

def geweke_diagnostic(chain, first=0.1, last=0.5):
    n = len(chain)
    first_end = int(n * first)
    last_start = int(n * (1 - last))
    mean_first = np.mean(chain[:first_end])
    mean_last = np.mean(chain[last_start:])
    var_first = np.var(chain[:first_end]) / first_end
    var_last = np.var(chain[last_start:]) / (n - last_start)
    z = (mean_first - mean_last) / np.sqrt(var_first + var_last)
    p_value = 2 * (1 - sp_stats.norm.cdf(abs(z)))
    return z, p_value

print("\n=== Geweke Diagnostic ===")
for name in ['beta', 'sigma']:
    samples = trace.posterior[name].values.flatten()
    z, p = geweke_diagnostic(samples)
    print(f"{name}: z={z:.3f}, p={p:.3f} {'✓' if abs(z) < 2 else '✗'}")

Related Topics

• See Bayesian Linear Regression for when MCMC is needed
• See Hierarchical Bayesian Models for complex models requiring MCMC
• See Bayesian Statistics for the Bayesian framework

Key Takeaways