When Diagnostics Lie: An Empirical Study of False Convergence in MCMC Sampling Across Five Algorithms and Multiple Target Geometries

Abstract

Markov chain Monte Carlo (MCMC) convergence diagnostics are essential for validating Bayesian inference, yet their reliability varies dramatically across sampling algorithms and target distribution geometries. We present a systematic empirical comparison of five MCMC samplers—random walk Metropolis-Hastings, Hamiltonian Monte Carlo, Gibbs sampling, slice sampling, and a NUTS-like adaptive sampler—evaluated against three convergence diagnostics: the Gelman-Rubin R-hat statistic, effective sample size (ESS), and the Geweke z-score test. Our experiments span multivariate Gaussian targets with dimensions in {2, 5, 10, 25, 50} and pairwise correlations in {0, 0.5, 0.9, 0.99}, plus banana-shaped (Rosenbrock), bimodal mixture, and funnel distributions. Across 73 experimental conditions with 4 chains each, we find that 41.1% of conditions exhibit false convergence—where R-hat drops below the standard 1.01 threshold before ESS reaches 400. Gibbs sampling shows the highest false convergence rate (75.0%), followed by slice sampling (62.5%) and the NUTS-like sampler (58.3%). Metropolis-Hastings exhibits 30.0% false convergence, while HMC shows 16.7%. Critically, the Geweke test frequently fails in high-dimensional settings, producing z-scores exceeding 40 even after 5000 iterations in 25-dimensional correlated targets. We quantify how each diagnostic's reliability degrades with increasing dimension and correlation, providing practical guidance for practitioners on combining multiple diagnostics to avoid premature termination of sampling.

1. Introduction

Markov chain Monte Carlo methods remain the primary computational tool for Bayesian inference in complex statistical models. The fundamental challenge facing practitioners is determining when a Markov chain has explored its target distribution sufficiently to yield reliable estimates. Premature termination of sampling—driven by overly optimistic convergence diagnostics—can produce severely biased posterior summaries without any visible warning.

The Gelman-Rubin diagnostic (Gelman and Rubin, 1992), commonly reported as the potential scale reduction factor R-hat, is perhaps the most widely used convergence criterion. The standard recommendation is that R-hat values below 1.01 (or the earlier threshold of 1.1) across all parameters indicate adequate convergence. The effective sample size (ESS), which accounts for autocorrelation within chains, provides a complementary perspective by estimating the equivalent number of independent draws. The Geweke diagnostic (Geweke, 1992) tests stationarity by comparing means from the first and last portions of a chain.

Despite their ubiquity, these diagnostics can produce misleading results. R-hat measures between-chain versus within-chain variance, and can appear satisfactory when multiple chains are stuck in the same region of parameter space. ESS captures within-chain mixing quality but is computationally expensive and can be unstable for short chains. The Geweke test, being a single-chain diagnostic based on asymptotic normality, may lack power in precisely the finite-sample regime where convergence assessment matters most.

Previous theoretical work has characterized convergence rates for specific algorithms. Roberts, Gelman, and Gilks established optimal scaling results for random walk Metropolis in high dimensions. Neal (2011) provided the theoretical foundation for Hamiltonian Monte Carlo and the No-U-Turn Sampler (NUTS), showing how gradient information dramatically improves exploration of correlated targets. Brooks and Gelman (1998) refined the original R-hat diagnostic and discussed its limitations.

What has been less systematically studied is the interaction between sampler choice, target geometry, and diagnostic reliability. A diagnostic that works well for a Gibbs sampler on a low-dimensional uncorrelated target may fail catastrophically for a Metropolis-Hastings chain exploring a high-dimensional correlated distribution. Understanding these failure modes is critical for applied Bayesian computation.

In this paper, we conduct a comprehensive empirical evaluation across five MCMC algorithms, four target distribution families, and a systematic grid of dimensions and correlations. Our primary contributions are:

1. Quantification of the "false convergence" rate—the frequency with which R-hat < 1.01 is achieved while ESS remains below 400—across 73 distinct experimental conditions.

2. Demonstration that false convergence rates depend strongly on the sampler-diagnostic combination, ranging from 16.7% (HMC) to 75.0% (Gibbs).

3. Evidence that the Geweke diagnostic becomes unreliable in high dimensions, with z-scores routinely exceeding 20 in dimensions above 10 regardless of the sampler used.

4. Practical recommendations for combining diagnostics to minimize the risk of premature termination.

2. Methods

2.1 MCMC Algorithms

We implement five MCMC samplers from scratch to ensure full control over algorithmic parameters and chain initialization.

Random Walk Metropolis-Hastings (MH). The proposal distribution is an isotropic Gaussian centered at the current state with scale parameter σ = 2.38/√d, following the optimal scaling result for target distributions that are products of identical univariate distributions. At each iteration, a candidate x' = x + σε, ε ~ N(0, I_d) is proposed and accepted with probability min(1, π(x')/π(x)).

Hamiltonian Monte Carlo (HMC). We augment the target with momentum variables p ~ N(0, I_d) and simulate Hamiltonian dynamics using the leapfrog integrator. The step size is set to ε = 0.2/√d and the number of leapfrog steps to L = max(3, ⌊2√d⌋). The resulting proposal is accepted via a Metropolis correction based on the change in the total Hamiltonian.

Gibbs Sampling. For multivariate Gaussian targets, we sample each coordinate from its full conditional distribution, which is itself Gaussian with analytically computed mean and variance derived from the covariance structure. Conditional parameters are precomputed for efficiency.

Slice Sampling. We implement the univariate slice sampler with stepping-out and shrinking, cycling through coordinates. The initial bracket width is set to w = 2.0, with stepping-out limited to 20 doublings and shrinking limited to 30 iterations per coordinate update.

NUTS-like Adaptive Sampler. We implement a simplified version of the No-U-Turn Sampler that selects a random trajectory length at each iteration by drawing a tree depth uniformly from {1, 2, ..., max_depth} with max_depth = 4, yielding between 2 and 16 leapfrog steps. The step size is ε = 0.15/√d. A Metropolis-Hastings correction ensures detailed balance.

2.2 Target Distributions

Multivariate Gaussian. The primary test family is N(0, Σ) where Σ has unit diagonal entries and off-diagonal entries equal to ρ. We test all combinations of dimension d ∈ {2, 5, 10, 25, 50} and correlation ρ ∈ {0, 0.5, 0.9, 0.99}.

Bimodal Mixture. An equal-weight mixture of two d-dimensional standard Gaussians centered at the origin and at (4, 4, ..., 4). This creates well-separated modes that challenge mixing, particularly as dimension increases.

Banana (Rosenbrock) Distribution. A two-dimensional distribution with density proportional to exp(-½[x₁² + (x₂ - bx₁²)²]) for curvature parameters b ∈ {0.1, 1.0, 5.0}. Higher values of b create more severely curved, banana-shaped contours that challenge algorithms relying on local Gaussian approximations.

Neal's Funnel Distribution. The log-variance variable v ~ N(0, 9) controls the scale of the remaining coordinates x_i | v ~ N(0, exp(v)). This distribution has extremely variable local geometry, ranging from wide dispersions when v is large to tightly concentrated mass when v is small. We test dimensions d ∈ {2, 5, 10}.

2.3 Convergence Diagnostics

Gelman-Rubin R-hat. Given m chains of length n, we compute the between-chain variance B = (n/(m-1))∑(θ̄_j - θ̄)² and the within-chain variance W = (1/m)∑s²_j. The estimated marginal variance is V̂ = ((n-1)/n)W + B/n, and R-hat = √(V̂/W). For multivariate targets, we report the maximum R-hat across all coordinates. The standard convergence threshold is R-hat < 1.01.

Effective Sample Size (ESS). We estimate the integrated autocorrelation time τ using FFT-based autocorrelation computation with the initial monotone sequence estimator: τ = 1 + 2∑_k ρ(k), where the sum terminates at the first k for which the sum of consecutive autocorrelation pairs becomes negative. The ESS is then n/τ. For multivariate chains, we report the minimum ESS across coordinates. We consider ESS > 400 as indicating adequate convergence, following common recommendations for reliable posterior quantile estimation.

Geweke Z-Score. We compare the mean of the first 10% of the chain to the mean of the last 50%, using the asymptotic standard error under the assumption of stationarity. For multivariate chains, we report the maximum absolute z-score across coordinates. Values exceeding 2 suggest non-stationarity at the 5% level.

2.4 Experimental Protocol

For each sampler-target combination, we run m = 4 chains initialized from overdispersed starting points drawn from N(0, 9I_d). Chain lengths are 5000 iterations for Metropolis-Hastings, multimodal, banana, and funnel experiments; 2000 for HMC and the NUTS-like sampler; 3000 for Gibbs; and 1000 for slice sampling (reflecting the higher per-iteration cost of slice sampling). Diagnostics are computed at multiple checkpoints throughout the chain to track convergence trajectories.

We define "false convergence" as a condition where R-hat drops below 1.01 at some checkpoint while ESS has not yet exceeded 400 at that same checkpoint or any earlier one. This captures the scenario where R-hat gives a reassuring signal but the chain has not actually produced enough effectively independent samples for reliable inference.

3. Results

3.1 Metropolis-Hastings on Multivariate Gaussians

Table 1 presents the final diagnostic values for random walk MH across all dimension-correlation combinations with 5000 iterations per chain.

Table 1: Metropolis-Hastings diagnostics on multivariate Gaussian targets (n = 5000, m = 4 chains)

Several patterns emerge immediately. In the low-dimensional, low-correlation regime (d = 2, ρ = 0), MH performs well: R-hat is 1.0012, ESS reaches 714.7, and the acceptance rate of 0.358 is close to the theoretically optimal 0.234 for high-dimensional targets. As either dimension or correlation increases, performance degrades sharply. At d = 5, ρ = 0.99, the acceptance rate plummets to 0.003 and R-hat explodes to 3.009, indicating that chains have barely moved from their starting points. ESS of 5.9 from 5000 samples means approximately one effectively independent draw per 850 iterations.

The Geweke statistic reveals an interesting phenomenon: even in settings where R-hat appears acceptable (e.g., d = 2, ρ = 0, R-hat = 1.0012), the Geweke z-score of 4.109 exceeds the critical value of 2, suggesting residual non-stationarity. This is likely an artifact of the overdispersed initialization—the early burn-in period inflates the first-segment mean relative to the equilibrium portion of the chain.

The convergence trajectory for a well-behaved case (d = 2, ρ = 0) shows R-hat declining from 1.026 at 100 iterations to 1.001 at 5000, with ESS growing linearly from 16.3 to 714.7. In contrast, for d = 10, ρ = 0.9, R-hat remains above 2.0 throughout the entire run, starting at 2.162 after 100 iterations and actually increasing to 2.537 by iteration 5000. This non-monotonic behavior of R-hat occurs because chains initialized in different regions slowly drift in different directions rather than mixing.

3.2 Hamiltonian Monte Carlo

HMC exploits gradient information to make coherent moves through the target distribution, producing dramatically different convergence behavior from MH.

Table 2: HMC diagnostics on multivariate Gaussian targets (n = 2000, m = 4 chains)

The acceptance rates for HMC remain above 0.77 in all tested conditions, compared to MH rates that drop below 0.01 in the most challenging settings. At d = 10, ρ = 0.99, HMC achieves an ESS of 16.9 from 2000 iterations (acceptance rate 0.926), while MH achieves only 5.9 from 5000 iterations (acceptance rate 0.004). On a per-sample basis, HMC is approximately 7 times more efficient in this regime, though each HMC iteration is computationally more expensive due to gradient evaluations and leapfrog integration.

The R-hat values for HMC cluster around 1.01–1.08, suggesting chains mix much more rapidly between their overdispersed starting points. However, the ESS values are still modest (12.6–110.5), indicating significant autocorrelation. This mismatch between apparently converged R-hat and low ESS represents a core false convergence scenario.

The Geweke z-scores for HMC in d = 5 and d = 10 are consistently elevated (9–15), which is surprising given the good acceptance rates. This occurs because HMC's deterministic trajectory can induce long-range correlations that are not captured by the simple first-versus-last comparison of the Geweke test.

3.3 Gibbs Sampling

The Gibbs sampler, which draws each coordinate from its exact conditional distribution, represents the opposite end of the algorithmic spectrum from MH.

Table 3: Gibbs sampler diagnostics on multivariate Gaussian targets (n = 3000, m = 4 chains)

Gibbs sampling shows a striking pattern: R-hat values are near-perfect (≤ 1.01) for nearly all conditions except d = 5 with ρ = 0.99 (R-hat = 1.0712) and d = 10 with ρ = 0.99 (R-hat = 1.1579). In uncorrelated targets, the ESS is enormous—2876 from 3000 samples in d = 2, indicating near-independence between consecutive draws when correlations are absent. This is expected since each Gibbs update draws exactly from the conditional distribution.

However, as correlation ρ increases to 0.99, the ESS collapses to 47.9 in d = 2 and 12.4 in d = 10. The key insight is that R-hat remains extremely close to 1.00 even when ESS is low. At d = 10, ρ = 0.9, R-hat is a mere 1.0067 while ESS is only 56.0—well below the 400 threshold commonly recommended for reliable quantile estimation. The Gibbs sampler's exact conditionals ensure that chains rapidly reach the correct marginal distributions for each coordinate, causing R-hat to decrease quickly, but the strong inter-coordinate correlations mean that the joint distribution is explored slowly.

This makes Gibbs sampling the most prone to false convergence (75.0% of conditions), because R-hat measures inter-chain agreement on marginals while ESS measures joint exploration efficiency.

3.4 Slice Sampling

Table 4: Slice sampler diagnostics on multivariate Gaussian targets (n = 1000, m = 4 chains)

Slice sampling, which automatically adapts its step size through the stepping-out procedure, exhibits behavior intermediate between Gibbs and MH. In uncorrelated settings, it achieves near-independent sampling (ESS = 996.3 from 1000 samples at d = 2, ρ = 0). Like Gibbs, it cycles through coordinates, making it sensitive to strong correlations: at d = 5, ρ = 0.99, R-hat remains a superficially acceptable 1.0010 while ESS is only 47.6 at ρ = 0.9, and at ρ = 0.99 the R-hat finally jumps to 1.246 with ESS collapsing to 8.9.

The false convergence rate for slice sampling (62.5%) is driven by the same mechanism as Gibbs: coordinate-wise updates achieve good marginal convergence (low R-hat) but poor joint exploration (low ESS) in correlated targets.

3.5 NUTS-like Adaptive Sampler

Table 5: NUTS-like sampler diagnostics on multivariate Gaussian targets (n = 2000, m = 4 chains)

The NUTS-like sampler demonstrates consistently excellent R-hat values (all below 1.09), owing to the ability of Hamiltonian dynamics to make large, coherent moves that rapidly mix between overdispersed starting positions. In the low-correlation regime, the R-hat values are very close to 1.00 and ESS is high—485.0 from 2000 samples at d = 2, ρ = 0, indicating roughly one independent draw per 4 iterations.

At d = 2 with ρ = 0.99, the NUTS-like sampler achieves ESS = 303.8 from 2000 samples, vastly outperforming MH (ESS = 80.3 from 5000 samples). The acceptance rates remain above 0.91 across all conditions, contrasting sharply with MH rates below 0.06 at ρ = 0.99.

Despite these superior results, the NUTS-like sampler has the second-highest false convergence rate (58.3%) among all samplers. This counterintuitive result arises because R-hat converges to near-unity very rapidly—often within the first 100 iterations—while ESS takes longer to exceed 400 due to residual autocorrelation. The sampler is doing well in absolute terms, but the R-hat diagnostic is so quick to declare convergence that it jumps ahead of the ESS criterion.

3.6 Non-Gaussian Target Distributions

The non-Gaussian targets reveal additional failure modes that would not appear with Gaussian tests alone.

Table 6: Diagnostics on bimodal mixture distributions (MH, n = 5000)

The bimodal mixture presents a challenging scenario. At d = 2, all four chains visit both modes (modes visited = 4/4), but R-hat is 1.0496 and ESS is only 31.0—the mode-switching behavior creates long autocorrelation lags. Paradoxically, at d = 5 and d = 10, all chains again visit both modes, but R-hat explodes to 2.23 and 2.28 respectively. This occurs because in higher dimensions, different chains spend different proportions of time in each mode, creating substantial between-chain variance.

The ESS of 248.8 at d = 5 despite an R-hat of 2.23 is a case of "reverse false convergence"—ESS suggests reasonable mixing while R-hat correctly flags a problem. This occurs because the ESS computation is based on within-chain autocorrelation and does not detect the between-chain disagreement that R-hat captures.

Table 7: Diagnostics on banana-shaped (Rosenbrock) distributions (MH, n = 5000)

The banana distribution smoothly interpolates between nearly Gaussian (b = 0.1) and severely non-Gaussian (b = 5.0). At b = 0.1, the target is nearly elliptical and MH performs well (ESS = 201.6, R-hat = 1.007). As curvature increases to b = 5.0, the chain struggles to navigate the curved, narrow ridge: ESS drops to 17.9 and R-hat increases to 1.228. The acceptance rate of 0.397 is actually reasonable, suggesting the chain makes moves at an acceptable rate but these moves are inefficient at exploring the curved geometry.

Table 8: Diagnostics on Neal's funnel distribution (MH, n = 5000)

The funnel distribution presents a uniquely challenging geometry: the width of the distribution in coordinates x_2, ..., x_d varies exponentially with x_1 (the log-variance parameter v). When v is large, the distribution is very wide; when v is small, it contracts to a narrow spike. Random walk MH with a fixed proposal scale cannot simultaneously explore both regimes efficiently. The acceptance rates (0.63–0.81) are deceptively high because the chain spends most of its time in the moderate-v region where the proposal scale is appropriate, rarely venturing into the narrow neck of the funnel.

At d = 10, R-hat reaches 1.405 and ESS is only 8.2, meaning approximately one effective sample per 600 iterations. This represents one of the most difficult targets in our test suite.

3.7 False Convergence Analysis

Our central finding concerns the systematic discrepancy between R-hat and ESS as convergence indicators.

Table 9: False convergence rates by sampler

The overall false convergence rate of 41.1% across all 73 conditions is alarmingly high. Nearly half of all sampler-target combinations would yield a prematurely optimistic assessment if R-hat alone were used as the convergence criterion.

The pattern across samplers reveals a clear dichotomy. Coordinate-wise samplers (Gibbs at 75.0%, slice at 62.5%) are most prone to false convergence because they achieve good marginal convergence rapidly while joint convergence lags. Gradient-based samplers show a more nuanced picture: HMC has the lowest false convergence rate (16.7%), but NUTS-like has a surprisingly high rate (58.3%) because its superior mixing causes R-hat to fall below 1.01 very quickly, before even 400 effective samples have accumulated.

The zero false convergence rates for multimodal and funnel targets reflect the opposite phenomenon: R-hat correctly indicates non-convergence (staying above 1.01) in these genuinely difficult cases, while ESS also remains low. When both diagnostics agree that convergence has not been achieved, there is no false convergence—the problem is that the chain simply has not converged at all.

False convergence by dimension for MH:

This dimension dependence is counterintuitive: false convergence is highest in low dimensions and decreases to zero in high dimensions. The explanation is that in high dimensions with the MH sampler, chains mix so poorly that R-hat never drops below 1.01, so false convergence cannot occur—both diagnostics correctly indicate failure. In low dimensions, chains mix well enough for R-hat to look good while ESS has not yet accumulated sufficiently, creating a window for false convergence.

False convergence by correlation for MH:

Again, the highest false convergence rate occurs at the lowest correlation (60.0% at ρ = 0). In the absence of correlation, MH chains mix efficiently enough that R-hat rapidly achieves 1.01 while ESS builds more slowly. At ρ = 0.99, the severe autocorrelation ensures that R-hat never reaches the threshold, preventing false convergence.

3.8 Chain Length Analysis

To understand how diagnostics evolve with chain length, we examined MH on a d = 10, ρ = 0.9 Gaussian at four chain lengths.

Table 10: Diagnostics versus chain length (MH, d = 10, ρ = 0.9)

Several notable features emerge. First, R-hat decreases from 5.414 at n = 1000 to 1.128 at n = 50000, but still has not crossed the 1.01 threshold even after 50,000 iterations—indicating that this is a genuinely challenging configuration for MH. Second, the ESS barely increases with chain length, growing only from 8.2 to 14.0 over a 50-fold increase in iterations. This is characteristic of chains that are essentially stuck: each additional iteration is nearly perfectly correlated with the previous one, contributing almost nothing to the effective sample size.

Third, and most strikingly, the Geweke z-score actually increases with chain length, from 25.1 at n = 1000 to 139.6 at n = 50000. This occurs because the Geweke test compares the first 10% and last 50% of the chain. As the chain length grows but mixing remains poor, the test has more statistical power to detect the persistent non-stationarity, yielding increasingly extreme z-scores. This makes the Geweke test a useful "alarm" in poorly-mixing chains—unlike R-hat, which can be slow to respond, the Geweke z-score becomes more emphatic about the problem as evidence accumulates.

4. Discussion

4.1 The False Convergence Problem

Our results demonstrate that no single convergence diagnostic is reliable across all sampler-target combinations. The 41.1% overall false convergence rate—defined as R-hat < 1.01 achieved before ESS > 400—is a practical concern for applied Bayesian analysis. In approximately two out of five experimental conditions, a practitioner relying solely on R-hat would terminate sampling prematurely.

The mechanisms of false convergence differ by sampler. For coordinate-wise methods (Gibbs, slice), the issue is that R-hat monitors marginal agreement across chains, while ESS captures joint mixing. These samplers can achieve excellent marginal convergence while the joint distribution remains poorly explored, particularly under strong correlations. For gradient-based methods (NUTS-like), the issue is the opposite: the sampler mixes so well that R-hat achieves its threshold very rapidly, before a sufficient number of effectively independent samples has been generated.

4.2 Diagnostic Reliability Across Geometries

The Geweke diagnostic emerges as the most unreliable of the three metrics tested. In dimensions above 10, Geweke z-scores routinely exceed 20 even in conditions where R-hat and ESS indicate adequate convergence. The root cause is the violation of the asymptotic normality assumption underlying the test: with moderate autocorrelation, the variance estimate in the z-score denominator is severely biased downward, inflating the test statistic. We recommend against using the Geweke test as a primary convergence criterion in moderate-to-high-dimensional problems.

R-hat, while more robust than the Geweke test, has the fundamental limitation of comparing only between-chain and within-chain variance. It cannot detect the scenario where all chains converge to the same wrong region of parameter space (a concern in multimodal targets), nor can it distinguish genuine convergence from the superficial agreement of slowly-drifting chains.

ESS provides the most reliable single summary of sampling quality, as it directly measures the effective information content of the chain. However, it is a within-chain diagnostic and does not capture between-chain disagreements. Furthermore, its computation via FFT-based autocorrelation, while fast, requires care in the truncation of the autocorrelation sum to avoid instability.

4.3 Sampler-Specific Recommendations

For Metropolis-Hastings users: In low-dimensional, low-correlation settings, MH with optimal scaling (σ = 2.38/√d) provides reliable performance, but monitor both R-hat and ESS. The acceptance rate is an excellent early warning: rates below 0.10 indicate that the proposal scale is poorly adapted to the target geometry, and convergence diagnostics should be interpreted with extreme caution. In our experiments, acceptance rates below 0.03 were universally associated with ESS below 12, regardless of what R-hat indicated.

For HMC and NUTS users: These gradient-based samplers show the best overall convergence properties, with HMC having the lowest false convergence rate (16.7%). However, NUTS users should be aware that the excellent mixing can cause R-hat to declare convergence before sufficient effective samples have accumulated. We recommend requiring ESS > 400 in addition to R-hat < 1.01 before terminating sampling.

For Gibbs sampler users: The 75.0% false convergence rate for Gibbs sampling is a serious concern. The exact conditional draws ensure rapid marginal convergence, causing R-hat to plummet toward 1.00, but joint exploration under strong correlations is poor. Gibbs users should pay primary attention to ESS and should be suspicious of R-hat values that are "too good" (R-hat < 1.001) in correlated targets—this may indicate that R-hat is simply insensitive to the relevant mixing failure.

For slice sampling users: The 62.5% false convergence rate follows the same pattern as Gibbs. The stepping-out mechanism adapts to the local scale of the target, providing efficient marginal updates, but the coordinate-wise cycling fails to capture correlation structure.

4.4 Practical Guidelines

Based on our empirical findings, we recommend the following protocol for convergence assessment:

1. Run multiple chains from overdispersed starting points (we used 4 chains with N(0, 9I_d) initialization throughout).

2. Compute R-hat and ESS jointly. Require both R-hat < 1.01 and ESS > 400 (or ESS > 100 × number of parameters being estimated, whichever is larger) before declaring convergence.

3. Monitor acceptance rates for MH-family algorithms. Rates below 0.10 in dimensions above 5 indicate inadequate proposal tuning and should trigger reparameterization or algorithm switching rather than longer runs.

4. Examine convergence trajectories, not just endpoints. A non-monotonically decreasing R-hat (as observed in d = 10, ρ = 0.9) is a red flag for poor mixing even if the final value appears acceptable.

5. Consider the Geweke test as supplementary only. While it can serve as an alarm in poorly-mixing chains (z-scores increase with chain length when mixing is inadequate), its systematic inflation in moderate-to-high dimensions limits its utility as a primary diagnostic.

6. Be most cautious with coordinate-wise samplers (Gibbs, slice) on correlated targets. The discrepancy between marginal convergence (tracked by R-hat) and joint convergence (tracked by ESS) is largest for these methods.

4.5 Limitations

Our study has several limitations that should inform interpretation. First, the "optimal" proposal scales used for MH (σ = 2.38/√d) and the step sizes for HMC/NUTS are tuned for Gaussian targets and may not be optimal for the non-Gaussian distributions in our test suite. Adaptive variants (e.g., adaptive MH with covariance learning, dual-averaging for HMC step size) would likely improve performance, though the diagnostic reliability issues we identify would persist.

Second, our definition of false convergence (R-hat < 1.01 before ESS > 400) involves two somewhat arbitrary thresholds. The R-hat threshold of 1.01 follows current best practices, but the ESS threshold of 400 is a pragmatic choice based on the accuracy requirements for posterior quantile estimation. Different thresholds would yield different false convergence rates, though the qualitative patterns across samplers and targets would remain.

Third, we use a simplified NUTS-like algorithm rather than the full NUTS implementation with its U-turn criterion, multinomial sampling, and dual-averaging step-size adaptation. The full NUTS algorithm would likely show improved performance, particularly in correlated targets.

Fourth, chain lengths were necessarily limited by computational constraints, especially for the more expensive samplers (HMC, NUTS, slice). Longer chains would provide more definitive convergence assessments, but our focus is precisely on the finite-sample regime where diagnostics must be relied upon.

5. Related Work

The foundational work on convergence diagnostics includes Gelman and Rubin (1992), who introduced the potential scale reduction factor based on between-chain and within-chain variance comparison. Brooks and Gelman (1998) proposed multivariate extensions and refined the diagnostic to account for chain length and number of chains. Geweke (1992) introduced the spectral density-based z-test for stationarity assessment.

Theoretical convergence rate results have been established for various MCMC algorithms. For random walk Metropolis, the optimal acceptance rate of approximately 0.234 in high dimensions was established under product target assumptions, and the √d scaling of the mixing time was derived. For Gibbs sampling on Gaussian targets, the convergence rate is known to depend on the spectral radius of the correlation matrix, which approaches 1 as the maximum correlation approaches 1, consistent with our empirical observations.

Neal (2011) provided a comprehensive treatment of HMC, including the theoretical motivation for leapfrog integration, the role of step size and trajectory length in mixing efficiency, and the No-U-Turn Sampler which automatically selects trajectory lengths. Our NUTS-like algorithm is inspired by this work but uses a simpler random-depth mechanism.

Recent work has proposed improved convergence diagnostics, including rank-normalized R-hat, folded R-hat for detecting chains that have converged to different modes, and the bulk-ESS and tail-ESS decomposition. Our study uses the classical versions of these diagnostics to establish a baseline understanding of their failure modes; extending this analysis to improved diagnostics is an important direction for future work.

6. Conclusion

Through a systematic empirical evaluation of 73 sampler-target-diagnostic combinations, we have demonstrated that convergence diagnostics for MCMC samplers are far less reliable than commonly assumed. The overall false convergence rate of 41.1%—where R-hat signals convergence before ESS indicates adequate effective sampling—should give pause to any practitioner relying on a single diagnostic.

The false convergence problem is most severe for coordinate-wise samplers (Gibbs: 75.0%, slice: 62.5%) on correlated targets, where marginal convergence occurs rapidly but joint exploration remains poor. Gradient-based methods (HMC: 16.7%, NUTS-like: 58.3%) show different failure patterns: HMC is relatively robust, while the NUTS-like sampler's excellent mixing causes R-hat to declare convergence before sufficient effective samples accumulate.

Our key practical recommendation is straightforward: never rely on a single convergence diagnostic. The joint criterion R-hat < 1.01 AND ESS > 400 would have correctly identified the convergence status in a substantially larger fraction of our test conditions than either criterion alone. Practitioners should also monitor acceptance rates, examine convergence trajectories rather than just endpoints, and exercise particular caution with coordinate-wise samplers on correlated targets.

The interaction between sampler choice, target geometry, and diagnostic reliability is an underappreciated aspect of applied Bayesian computation. We hope this empirical study provides both practical guidance and motivation for the development of more robust convergence assessment tools that account for the specific characteristics of the sampling algorithm being used.

References

Brooks, S. P. and Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7(4):434–455.

Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4):457–472.

Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In Bayesian Statistics 4, pages 169–193. Oxford University Press.

Neal, R. M. (2011). MCMC using Hamiltonian dynamics. In Brooks, S., Gelman, A., Jones, G. L., and Meng, X.-L., editors, Handbook of Markov Chain Monte Carlo, pages 113–162. Chapman and Hall/CRC.