PaperlayerView exact revision on arXiv

Source-linked AI summary

Convergence diagnostics for Markov chain Monte Carlo

Vivekananda Roy

arXiv:1909.11827v2stat.CO

TL;DR

MCMC users must determine suitable starting points and stopping times while empirical convergence diagnostics may provide unreliable evidence. This review surveys diagnostic tools and theoretically grounded stopping rules, concluding that fixed-width and effective-sample-size rules are preferable to some empirical alternatives.

  • Problem

    MCMC implementation requires determining where to start and when to stop so that chains and Monte Carlo estimators have sufficiently converged.

  • Method

    The review surveys convergence diagnostics, theoretically grounded stopping rules, and their behavior in three illustrative MCMC examples.

  • Results

    Some empirical diagnostics may prematurely terminate simulations, whereas the review recommends fixed-width and effective-sample-size-based stopping rules.

  • Takeaways & Limitations

    Practitioners should avoid relying purely on empirical diagnostics, particularly when the target distribution may have multiple modes.

  • Takeaways & Limitations

    Empirical diagnostics can falsely indicate convergence when chains have not run long enough to move between modes.

Abstract

from arXiv · show

Markov chain Monte Carlo (MCMC) is one of the most useful approaches to scientific computing because of its flexible construction, ease of use and generality. Indeed, MCMC is indispensable for performing Bayesian analysis. Two critical questions that MCMC practitioners need to address are where to start and when to stop the simulation. Although a great amount of research has gone into establishing convergence criteria and stopping rules with sound theoretical foundation, in practice, MCMC users often decide convergence by applying empirical diagnostic tools. This review article discusses the most widely used MCMC convergence diagnostic tools. Some recently proposed stopping rules with firm theoretical footing are also presented. The convergence diagnostics and stopping rules are illustrated using three detailed examples.

1. INTRODUCTION

MCMC implementation requires deciding where to start collecting samples and when to stop, corresponding to convergence of the chain and Monte Carlo estimators. Although theoretical bounds and sample-size calculations can guide these decisions, finite simulations require assessing estimator accuracy without wasting resources.

  • Starting and stopping: MCMC practitioners must determine where to start and when to stop the algorithm, addressing chain convergence to stationarity and estimator convergence to population quantities.Under standard conditions, the chain distribution converges to its stationary distribution for any initial value as n →∞.
  • Monte Carlo estimation: Target-density quantities can often be written as expectations Eπg, which MCMC estimates with the post-burn-in average ¯g_n′,n.The estimator averages g(X_i) from i=n′+1 through n over n−n′ samples; setting n′=0 means no burn-in.
  • Stopping criteria: Finite simulations should stop only when the estimator ¯g_n′,n* has sufficiently converged to Eπg, because premature termination can produce inaccurate inference.The quality of the time-average estimator depends on the quality of the samples.
  • Stopping criteria: Unnecessarily long chains are also undesirable because they consume computational resources.Thus, convergence assessment must balance avoiding premature termination with limiting wasted computation.
  • Theoretical guidance: Theoretical analysis can provide rigorous stopping guidance through upper bounds on distance to stationarity and sample-size calculations based on asymptotic Monte Carlo error.These approaches can yield a burn-in value n′ or an honest stopping value n* when the necessary theoretical analysis is available.

2. MCMC diagnostics

This section describes diagnostic tools for deciding whether Markov chains have reached stationarity and stopping rules for ending MCMC sampling while using resources prudently.

  • 2. MCMC diagnostics: MCMC diagnostics assess convergence to stationarity and help determine when to stop sampling.Although longer chains generally improve Monte Carlo estimates, stopping rules support prudent resource use.

2.1. Honest MCMC

Honest MCMC defines rigorous burn-in and stopping procedures using total-variation convergence and Markov-chain central limit theory. These procedures rely on estimating Monte Carlo uncertainty, while drift-and-minorization techniques and related variance estimators support practical implementation.

  • Burn-in: The honest burn-in is the smallest iteration n′ meeting a predetermined total-variation precision threshold for convergence to π.The commonly illustrated cutoff 0.01 is arbitrary; any predetermined precision level may be used.
  • Burn-in: Total-variation convergence is generally unavailable, and constructing quantitative bounds is often difficult, although drift and minorization provides a key tool.Drift and minorization has been used to analyze a variety of MCMC algorithms.
  • Stopping rules: A Markov-chain central limit theorem makes the sample average a reliable estimator with standard error bσg,n/√n and supports confidence intervals for Eπg.The standard error should be reported with the point estimate because it indicates estimate reliability.
  • Stopping rules: The honest stopping rule requires a Markov-chain CLT and a consistent estimator of σg, with geometric ergodicity commonly used to establish both.The CLT requires total-variation convergence at a suitable rate.
  • Extensions: Stopping rules extend to vector-valued functions and quantile estimation, and sequential checks can be performed every l iterations to reduce computational burden.For vector-valued functions, a consistent covariance estimator yields an asymptotic confidence region.

2.2. Relative fixed-width stopping rules

Relative fixed-width stopping rules terminate MCMC simulations when estimated error is small relative to the estimated quantity or its standard deviation. In Bayesian and high-dimensional settings, relative SDFWSR is advocated, while multivariate settings motivate alternatives to rules based on marginal chains.

  • Relative fixed-width stopping rules: Relative magnitude FWSR stops after ñ iterations when the error bound is at most ε times the estimator ḡ_n.The rule uses t* bσ_g,n n^-1/2 + n^-1 ≤ ε ḡ_n.
  • Relative fixed-width stopping rules: Relative SDFWSR stops after ñ iterations when the error bound is at most ε times the estimated standard deviation.Its criterion replaces the estimator with bλ_g,n on the right-hand side.
  • Relative fixed-width stopping rules: In Bayesian applications, Flegal and Gong advocate relative SDFWSR, which Gong and Flegal prefer over FWSR based on marginal chains for high-dimensional functions without known Eπg magnitude.Here g is R^p-valued and p is large.
  • Relative fixed-width stopping rules: For multivariate settings, stopping rules based on p marginal chains may be inappropriate, motivating a proposed alternative stopping rule.The supplied passage indicates that the proposal stops the simulation but does not include its full criterion.

2.3. Effective sample size

Effective sample size (ESS) measures the number of independent draws equivalent to correlated Markov chain samples in estimator precision. Its estimated univariate or multivariate form can guide termination decisions and compare MCMC algorithms by computational and statistical efficiency.

  • Definition: ESS is the number of independent samples having the same standard error as a set of correlated Markov chain samples.The common definition is based on the chain length and the summed autocorrelations of g(X0) and g(Xi).
  • Multivariate ESS: For Rp-valued functions g, multivariate ESS (mESS) extends ESS to the multivariate setting using the population covariance matrix Λg.Vats et al. (2019) define mESS for multivariate output.
  • Stopping rule: MCMC simulation can terminate when a consistent estimator of ESS or mESS exceeds a prespecified threshold.The ESS stopping rule can target a desired precision for the volume of an asymptotic confidence region.
  • Efficiency comparison: Estimated ESS or mESS per unit time compares MCMC algorithms with the same stationary distribution in computational and statistical efficiency.ESS is implemented in R packages including coda and mcmcse.

2.4. Gelman-Rubin diagnostic

The Gelman–Rubin diagnostic assesses MCMC convergence by comparing variance estimates across multiple over-dispersed chains. Simulation is typically stopped when the potential scale reduction factor is sufficiently close to one, with 1.1 generally used as the cutoff; extensions address single-chain and multivariate settings.

  • Multiple-chain diagnostic: The GR diagnostic is widely used and relies on multiple chains initialized from a distribution over-dispersed relative to the target density.In practice, initial points are often chosen ad hoc.
  • Variance comparison: The diagnostic estimates within-chain and between-chain variance, combines them into a pooled variance estimate, and compares their ratio through the potential scale reduction factor.The between-chain estimate uses differences between chain means and the overall mean.
  • Stopping rule: 1.1 is the generally used cutoff for stopping simulation when the potential scale reduction factor is sufficiently close to one.Over-dispersed initialization makes the numerator tend to overestimate the target variance and the denominator tend to underestimate it, producing ˆR > 1 in finite samples.
  • Modified diagnostic: The modified GR statistic connects the diagnostic with effective sample size and permits computation from a single chain.Its expression differs slightly from the original Gelman–Rubin definition, although the widely used version remains common in practice.
  • Multivariate extension: The multivariate PSRF extends the diagnostic to multivariate convergence, using pooled, within-chain, and between-chain covariance matrices and stopping when ˆRp ≈1.Peltonen et al. propose visualization based on linear and discriminant component analysis to complement these diagnostics.

2.5. Two spectral density-based methods

This section presents Geweke’s and Heidelberger–Welch’s convergence diagnostics, which use spectral density estimates or related asymptotic variance estimates to assess stationarity. Both are univariate tools implemented in the CODA package.

  • Geweke diagnostic: Geweke’s diagnostic uses the spectral density at zero, Sg(0), as the asymptotic variance for estimating Eπg.It compares means from two parts of the chain while accounting for sample autocorrelation in the standard error.
  • Geweke diagnostic: Geweke suggests setting nA = 0.1n and nB = 0.5n for the two chain segments.The resulting Z score assumes independence between the segments and tests equality of their means.
  • Heidelberger–Welch diagnostic: Heidelberger–Welch’s diagnostic uses spectral density estimates and a Cramer-von Mises stationarity test based on an asymptotic Brownian-bridge assumption.The test statistic is the integral of Bn(t)^2 over t from 0 to 1.
  • Heidelberger–Welch diagnostic: The Heidelberger–Welch stationarity test is applied successively after rejecting the first 10%, 20%, and later portions of samples, up to 50%.Testing stops when the stationarity test is passed or 50% of the samples have been rejected.
  • Implementation: Both spectral density-based tools are implemented in the CODA package and are univariate diagnostics.The passage identifies CODA as the implementation package for both methods.

2.6. Raftery-Lewis diagnostic

The Raftery–Lewis diagnostic estimates burn-in and run length for quantile estimation, targeting probability accuracy ϵ with probability 1−α. It analyzes a binary indicator process and should be repeated across different quantiles because the diagnostic depends on q.

  • 2.6. Raftery-Lewis diagnostic: The method estimates burn-in and selects a run length so the probability estimate lies within [q − ϵ, q + ϵ] with probability 1 − α.It targets estimating u such that Pπ(g(X) ≤ u) = q.
  • 2.6. Raftery-Lewis diagnostic: The diagnostic transforms the target into the binary process W_n = I(g(X_n) ≤ u) and analyzes its two-state empirical transition matrix.Burn-in is estimated using eigenvalue analysis of the empirical transition matrix.
  • 2.6. Raftery-Lewis diagnostic: The diagnostic should be repeated for different quantiles because its results depend on q.The prescribed initial size is based on the standard asymptotic sample-size calculation for a Bernoulli(q) population.

2.7. Kernel density-based methods

Kernel density-based diagnostics assess convergence by comparing estimated whole distributions rather than summary moments. The methods include distance-based diagnostics between chains, KL-based tools for multimodal targets, and visualization for investigating chain clustering.

  • 2.7. Kernel density-based methods: Kernel density diagnostics declare convergence when the distance between density estimates from two chains or chain segments is close to zero.They assess the convergence of whole distributions, unlike the GR diagnostic, which compares summary moments.
  • 2.7. Kernel density-based methods: L1, Hellinger, and symmetric KL distances provide kernel density-based convergence diagnostics, with the KL tools introduced by Dixit and Roy (2017).The KL estimate is the average of the two directional divergences between density estimates.
  • 2.7. Kernel density-based methods: For multiple chains, Tool 1 uses the maximum estimated symmetric KL divergence across chain pairs and applies hypothesis-testing cutoffs.For multivariate chains, Dixit and Roy recommend marginal assessment with Bonferroni-adjusted significance levels when high-dimensional density estimation is unreliable.
  • 2.7. Kernel density-based methods: Tool 2 compares the chain’s kernel density estimate with the target density to detect divergence when chains may be stuck at the same mode.If T*2 > 0.05, the chain is interpreted as not having captured the target distribution adequately; numerical integration prevents use in high dimensions.
  • 2.7. Kernel density-based methods: A tile-plot visualization complements multi-chain diagnostics by marking chain pairs as same-cluster or different-cluster according to the KL cutoff.Applied marginally to multivariate chains, it helps identify variables responsible for inadequate mixing.

2.8. Graphical methods

Graphical MCMC diagnostics use trace, autocorrelation, and running mean plots to assess mixing, dependence, convergence, and whether simulation can stop. In multivariate settings, marginal plots do not reveal correlations among components, so cross-variable correlation must also be investigated.

  • Trace plots: Trace plots display each iteration’s chain realization against iteration number to visualize movement through the state space and assess mixing.They are used to diagnose convergence to stationarity.
  • Autocorrelation plots: Autocorrelation plots show lag-k sample autocorrelations; rapidly declining values indicate fast mixing, whereas persistently high values indicate strong dependence.The lag-k autocorrelation measures correlation between samples k steps apart.
  • Running mean plots: Running mean plots display Monte Carlo time-average estimates across iterations and should stabilize before the simulation is stopped.A non-converging running mean indicates that the simulation cannot yet be stopped.
  • Multivariate limitation: In multivariate analyses, marginal trace, autocorrelation, and running mean plots omit correlations among components, requiring separate checks for high cross-correlation.The individual plots are generally based on realizations of each marginal chain.

3. Examples

Three examples show that MCMC convergence diagnostics can indicate convergence prematurely or falsely, especially when chains mix slowly or remain trapped in a local mode. A real-data analysis illustrates how burn-in, PSRF, KL divergence, confidence-interval half-widths, and mESS can be used together to assess convergence and Monte Carlo accuracy.

  • Exponential target: For the exponential target, the θ = 5 chain mixes slowly, as shown by frequent flat trace segments and high autocorrelation.The θ = 0.5 and θ = 5 independence Metropolis chains illustrate burn-in and stopping-time choices alongside other diagnostics.
  • Exponential target: The θ = 5 chain’s final estimate is 0.778 rather than the truth EπX = 1, even after 300,000 iterations.Because a Markov chain CLT is unavailable for θ > 1, asymptotic confidence intervals cannot be computed for this chain.
  • Exponential target: PSRF falls below 1.1 before 100 iterations for θ = 0.5, causing premature termination despite the diagnostic’s apparent convergence.This example demonstrates that empirical diagnostics can falsely indicate convergence to stationarity or convergence of Monte Carlo estimates.
  • Sixmodal target: In the sixmodal example, PSRFs and MPSRF fall below 1.1 while all four chains remain stuck at the same local mode, falsely detecting convergence.Adaptive bivariate density plots also fail to reveal non-convergence because the chains look similar.
  • Real-data analysis: For the Anguilla logistic model, mESS is 55,775 versus a cutoff of 55,191 for ε = 0.02, but the ε = 0.01 cutoff is 220,766.The analysis removes 70,000 burn-in iterations, then runs each of three chains for 15,000 additional iterations; MPSRF is 1.004.

4. Conclusions and discussion

The article distinguishes diagnostics for convergence to stationarity, which inform burn-in, from diagnostics for sample-average convergence, which inform simulation termination. It cautions that empirical tools can fail or terminate prematurely, recommending fixed-width and ESS-based rules alongside careful handling of thinning, multimodality, and target propriety.

  • Conclusions and discussion: Diagnostics for stationarity convergence help determine burn-in, whereas diagnostics for sample-average convergence guide termination of the MCMC simulation.Analytical total-variation bounds for honest burn-in can be difficult to obtain or overly conservative.
  • Conclusions and discussion: Empirical diagnostics may prematurely terminate simulations and produce inferences far from the truth; fixed-width and ESS-based stopping rules are recommended.Many quantitative convergence diagnostics assume a Markov chain central limit theorem.
  • Conclusions and discussion: Thinning should be used only when storage is limited or evaluating functions is more expensive than sampling, because it wastes samples.Convergence diagnostics can be applied to thinned samples if thinning is used.
  • Conclusions and discussion: With multiple modes, parallel-chain or split-chain diagnostics can miss non-convergence when chains do not explore distinct high-density regions or move between modes.The authors recommend single long runs for final inference, noting that longer runs may reveal new support regions.
  • Conclusions and discussion: Practitioners should not rely purely on empirical diagnostics because they cannot certify convergence and may fail to flag an improperly assumed proper target density.This caution is especially important when multiple modes are suspected.
Loading 1909.11827v2…