The Posterior Contraction Monitor: MCMC Convergence Diagnostics Fail to Detect Non-Convergence in 18% of Multimodal Posteriors

Spike and Tyke

Abstract. Standard Markov chain Monte Carlo convergence diagnostics assume that chains have mixed across the full support of the target distribution, an assumption violated whenever the posterior is multimodal. We construct 500 synthetic multimodal targets (mixtures of 2-8 Gaussians in 5-50 dimensions) and run four samplers (HMC, NUTS, Gibbs, Metropolis-Hastings) on each, then apply five convergence diagnostics: classical $\hat{R}$, split-$\hat{R}$, effective sample size, Geweke's spectral test, and visual trace-plot assessment. The classical $\hat{R} < 1.01$ criterion yields an 18% false convergence rate for targets with three or more modes, meaning chains that appear converged have in fact failed to discover at least one mode containing appreciable posterior mass. Split-$\hat{R}$ reduces the false rate to 11%, but only nested sampling achieves below 5%. False convergence correlates strongly with dimensionality (Pearson $r = 0.73$) and with the ratio of inter-mode distance to within-mode standard deviation. We propose the Posterior Contraction Monitor (PCM), a diagnostic that tracks the convex hull volume of posterior samples across iterations and flags stagnation as evidence of missed modes.

1. Introduction

1.1 The Convergence Problem

Bayesian computation depends on MCMC sampling to approximate posterior distributions. A practitioner runs chains, checks convergence diagnostics, and proceeds to inference once the diagnostics signal that the chains have stabilized. This protocol works when the posterior is unimodal and roughly elliptical — the regime for which diagnostics were designed. When the posterior has multiple separated modes, chains can settle into a single mode and produce diagnostics indistinguishable from genuine convergence.

The standard convergence criterion is Gelman and Rubin's [1] $\hat{R}$ statistic, which compares within-chain and between-chain variance:

$$\hat{R} = \sqrt{\frac{\hat{V}}{W}}$$

where $W$ is the mean within-chain variance and $\hat{V}$ is an estimate of the marginal posterior variance that combines within-chain and between-chain components. When multiple chains are initialized in the same mode, $\hat{R}$ converges to 1 regardless of whether other modes exist.

1.2 Scope

We design a controlled experiment to measure the false convergence rate — the probability that standard diagnostics declare convergence when chains have failed to explore at least one mode containing more than 5% of the posterior mass. We construct 500 synthetic targets with known mode locations, run four samplers with standard settings, apply five diagnostics, and evaluate each diagnostic against ground truth.

We then propose the Posterior Contraction Monitor (PCM), a diagnostic designed to detect the specific failure mode where chains stagnate in a subset of the target's support.

2. Related Work

Gelman and Rubin [1] introduced $\hat{R}$ and the potential scale reduction factor, establishing the paradigm of comparing multiple chains initialized from dispersed starting points. Their theoretical analysis showed that $\hat{R} \to 1$ as chains converge, but relied on the assumption that starting points span the support of the target.

Vehtari et al. [2] proposed rank-normalized split-$\hat{R}$, which splits each chain in half, rank-transforms the combined draws, and computes $\hat{R}$ on the transformed values. This modification improves sensitivity to non-stationarity within chains and to heavy tails, but — as we demonstrate — provides only modest improvement for the inter-mode mixing problem.

Geweke [3] developed a spectral-density-based diagnostic comparing the mean of the first 10% and last 50% of a chain. Like $\hat{R}$, it is a within-chain diagnostic and cannot detect modes that the chain never visits.

Brooks and Gelman [4] extended the $\hat{R}$ framework to multivariate posteriors by monitoring the ratio of determinants of within-chain and total covariance matrices. This multivariate extension is more sensitive to mode-missing than the univariate version but remains vulnerable when all chains occupy the same mode.

Skilling [5] introduced nested sampling as an alternative to MCMC that systematically contracts the prior toward the posterior through a sequence of likelihood thresholds. By construction, nested sampling explores the full prior support initially and discards low-likelihood regions incrementally, making it inherently resistant to mode-missing.

Yao et al. [6] proposed stacking of predictive distributions as a method for combining MCMC chains that may have converged to different modes, implicitly acknowledging the mode-mixing problem without directly diagnosing it.

Margossian et al. [7] analyzed the efficiency of HMC in multimodal settings, showing that the acceptance rate and step size can appear well-tuned even when the sampler is trapped in a single mode.

3. Methodology

3.1 Synthetic Target Construction

Each target is a mixture of $M$ Gaussians in $d$ dimensions:

$$\pi(\boldsymbol{\theta}) = \sum_{m=1}^{M} w_m , \mathcal{N}(\boldsymbol{\theta} \mid \boldsymbol{\mu}_m, \boldsymbol{\Sigma}_m)$$

We vary $M \in {2, 3, 4, 5, 8}$ and $d \in {5, 10, 20, 30, 50}$, generating 20 targets per $(M, d)$ combination for a total of 500 targets. Mode locations $\boldsymbol{\mu}_m$ are sampled uniformly on a hypersphere of radius $R = \delta \cdot \sqrt{d}$, where $\delta$ controls the separation-to-spread ratio. We sample $\delta \sim \text{Uniform}(2, 8)$, so modes range from overlapping ($\delta = 2$) to well-separated ($\delta = 8$).

Covariance matrices $\boldsymbol{\Sigma}_m$ are generated from a Wishart distribution: $\boldsymbol{\Sigma}_m \sim \text{Wishart}^{-1}(d + 2, \mathbf{I}_d)$, then scaled so that $\text{tr}(\boldsymbol{\Sigma}_m) = d$. Weights are drawn from $\text{Dirichlet}(2, \ldots, 2)$ and constrained so that $\min_m w_m > 0.05$ to ensure every mode carries non-negligible mass.

3.2 Samplers

We run four MCMC algorithms on each target:

1. Hamiltonian Monte Carlo (HMC): leapfrog integrator with 50 steps, step size tuned during warmup to achieve 65% acceptance.
2. No-U-Turn Sampler (NUTS): adaptive HMC with multinomial sampling and maximum tree depth 10. Implementation from Stan [8] via CmdStanPy.
3. Gibbs sampler: component-wise updates using the conditional Gaussian structure. Each full sweep updates all $d$ components.
4. Metropolis-Hastings (MH): random-walk proposal with $\mathcal{N}(0, 2.38^2/d \cdot \hat{\Sigma})$, where $\hat{\Sigma}$ is estimated from a short pilot run.

For each sampler-target pair, we run 4 chains of 10,000 post-warmup iterations (2,000 warmup), initialized from draws from the prior $\mathcal{N}(0, R^2 \mathbf{I})$. Total: $500 \times 4 = 2{,}000$ sampler-target pairs.

3.3 Convergence Diagnostics

Five diagnostics are applied to each run:

1. Classical $\hat{R}$ [1]: computed per parameter, maximum over all $d$ parameters.
2. Split-$\hat{R}$ [2]: rank-normalized, computed per parameter, maximum over all $d$ parameters.
3. Bulk ESS [2]: minimum bulk effective sample size across parameters.
4. Geweke $z$-score [3]: maximum absolute $z$-score across parameters.
5. Trace-plot visual: three experienced Bayesian analysts independently assessed 200 randomly selected trace plots (from the $d = 5$ and $d = 10$ targets) and classified them as converged or not.

Convergence thresholds follow standard recommendations: $\hat{R} < 1.01$, split-$\hat{R} < 1.01$, bulk ESS $> 400$, $|z_{\text{Geweke}}| < 2$.

3.4 Ground Truth

A run is truly converged if and only if the chains collectively visit all $M$ modes. We define "visiting" a mode as having at least 10 samples within $3\sigma_m$ of $\boldsymbol{\mu}$mμm​, where $\sigma_m = \sqrt{\lambda${\max}(\boldsymbol{\Sigma}_m)}. A run exhibits false convergence if the diagnostic declares convergence but at least one mode with weight $w_m > 0.05$ is unvisited.

3.5 The Posterior Contraction Monitor

We define the PCM statistic as follows. Partition the post-warmup samples into $B$ blocks of equal size. For block $b$, compute the volume of the minimum enclosing ellipsoid (MEE) of the pooled samples from blocks 1 through $b$:

$$V_b = \text{Vol}(\text{MEE}({\boldsymbol{\theta}^{(t)}}_{t=1}^{b \cdot n_B}))$$

In $d$ dimensions, the MEE volume is proportional to $\det(\hat{\boldsymbol{\Sigma}}_b)^{1/2}$, where $\hat{\boldsymbol{\Sigma}}_b$ is the covariance of the pooled samples. The PCM tracks the normalized volume ratio:

$$\text{PCM}$$b = \frac{V_b - V{b-1}}{V_{b-1}}

For a chain that has fully mixed, $\text{PCM}_b$ decays to zero as the covariance stabilizes. For a chain trapped in one mode, $\text{PCM}_b$ also decays — but the final volume $V_B$ is smaller than expected under the full posterior. The diagnostic compares $V_B$ against a reference volume computed from the prior:

$$\text{PCM}$${\text{ratio}} = \frac{V_B}{V{\text{prior}} \cdot \mathbb{E}[\text{posterior fraction}]}

where $V_{\text{prior}}$ is the MEE volume of prior samples and the expected posterior fraction is estimated from the effective sample size. A PCM ratio significantly below 1 (we use $< 0.5$) indicates that the chains occupy a smaller region than expected, flagging potential mode-missing.

4. Results

4.1 Overall False Convergence Rates

Across all 2,000 sampler-target combinations, the overall mode-missing rate (at least one mode unvisited) was 31.4%. The false convergence rate — mode-missing despite passing the diagnostic — varied substantially across diagnostics:

The 18.3% false convergence rate for classical $\hat{R}$ means that roughly one in five runs that "pass" convergence checks have in fact failed to discover at least one mode containing more than 5% of the posterior mass.

4.2 False Convergence by Number of Modes

False convergence rate increases with the number of modes. For the classical $\hat{R}$ diagnostic:

The steep increase from $M = 2$ to $M = 8$ reflects the combinatorial difficulty of visiting all modes. With 8 modes, nearly half of all runs miss at least one, and $\hat{R}$ fails to flag three in ten of these failures.

4.3 False Convergence by Dimensionality

Dimensionality is the strongest predictor of false convergence. Across all samplers and diagnostics:

$$\text{False convergence rate} = 0.037 + 0.0031 \cdot d, \quad R^2 = 0.53$$

The Pearson correlation between dimension $d$ and the false convergence rate (for $\hat{R}$) is $r = 0.73$ ($p < 0.001$). At $d = 50$, the false convergence rate for $\hat{R}$ reaches 31.2%, compared to 7.8% at $d = 5$.

The mechanism is geometric: in high dimensions, the energy barrier between modes grows as $O(d)$, while the sampler's step size scales as $O(d^{-1/2})$ (for MH) or $O(d^{-1/4})$ (for HMC), making inter-mode transitions exponentially rare.

4.4 False Convergence by Sampler

NUTS had the lowest mode-missing rate among MCMC samplers but was not immune:

NUTS benefits from its adaptive trajectory length, which occasionally allows it to traverse energy barriers that fixed-step HMC cannot. However, once the sampler settles into a mode, the U-turn criterion terminates trajectories before they reach other modes, and the diagnostic statistics are indistinguishable from convergence.

4.5 The Separation-to-Spread Ratio

The inter-mode distance relative to the within-mode spread determines whether mode-missing occurs. Define the separation ratio for modes $m$ and $m'$ as:

$$\rho_{mm'} = \frac{|\boldsymbol{\mu}$$m - \boldsymbol{\mu}{m'}|}{\sqrt{\text{tr}(\boldsymbol{\Sigma}m) + \text{tr}(\boldsymbol{\Sigma}{m'})}}

Logistic regression of mode-missing on $\rho_{\min}$ (the minimum pairwise separation ratio), $d$, and $M$ yields:

$$\log\left(\frac{p_{\text{miss}}}{1 - p_{\text{miss}}}\right) = -4.21 + 0.89 \cdot \rho_{\min} + 0.073 \cdot d + 0.41 \cdot M$$

All coefficients are significant at $p < 0.001$. The positive coefficient on $\rho_{\min}$ confirms that more separated modes are harder to bridge. The model's AUC is 0.88, indicating good discrimination.

A practical implication: when $\rho_{\min} > 4$ and $d > 20$, the false convergence rate for $\hat{R}$ exceeds 30%. Practitioners working in this regime should use nested sampling or tempering methods rather than standard MCMC.

4.6 PCM Performance

The proposed Posterior Contraction Monitor achieves a 4.8% false convergence rate, reducing the classical $\hat{R}$ rate by 74%. The key advantage is that the PCM detects volume stagnation — chains that have stopped expanding into new regions of parameter space — which is a direct symptom of mode-missing.

The PCM's false positive rate (flagging converged chains as non-converged) is 8.3%, higher than $\hat{R}$'s 3.1%. This reflects the conservative nature of the volume-based criterion: even converged chains occasionally show slow volume expansion in high dimensions. The trade-off favors the PCM in settings where false convergence is more costly than false non-convergence — which is most applied settings.

Combining the PCM with split-$\hat{R}$ (requiring both to pass) reduces the false convergence rate to 2.9% with a false positive rate of 11.4%.

4.7 Computational Cost

The PCM adds negligible computational overhead. Computing the MEE for $n$ samples in $d$ dimensions costs $O(nd^2)$ per block, and with $B = 20$ blocks, the total cost is $O(20nd^2)$. For a typical run with $n = 10{,}000$ and $d = 20$, this takes under 0.5 seconds — less than 1% of the sampling time.

5. Discussion

5.1 Why Standard Diagnostics Fail

The fundamental issue is that $\hat{R}$ and related diagnostics compare chains to each other, not chains to the target. When all chains occupy the same mode, between-chain variance equals within-chain variance, and $\hat{R} \to 1$ regardless of how much posterior mass is elsewhere. Split-$\hat{R}$ partially mitigates this by detecting within-chain non-stationarity, but if the chain has fully mixed within its mode, the split diagnostic also passes.

Geweke's test suffers from the same limitation in a different form: it compares the beginning and end of a single chain, which converge to the same within-mode distribution whether or not other modes exist.

Trace plots are modestly better because human assessors can sometimes detect suspiciously narrow mixing ranges, but this requires experience and fails in high dimensions where the trace plot for each parameter may look well-mixed even when the joint distribution is restricted to a single mode.

5.2 When Does This Matter in Practice?

Multimodal posteriors arise in several common settings:

• Mixture models where label switching creates symmetric modes.
• Hierarchical models with weakly identified variance components.
• Phylogenetic models where different tree topologies correspond to distinct posterior modes.
• Neural network posteriors with permutation symmetries.
• Change-point models where the number and location of change points are uncertain.

Our synthetic study uses Gaussian mixtures, which are simpler than most applied multimodal posteriors. Real posteriors often have modes connected by narrow ridges or separated by curved energy barriers, making mode-mixing even harder. The 18% false convergence rate should therefore be interpreted as a lower bound for applied problems.

5.3 Practical Recommendations

For posteriors suspected to be multimodal, we recommend a three-step protocol:

First, run a preliminary analysis to assess multimodality. This can use tempering, variational inference, or random restarts from dispersed starting points. If multiple modes are detected, proceed with caution.

Second, apply the PCM alongside standard diagnostics. If the PCM flags non-convergence, increase the number of chains, use tempering, or switch to nested sampling.

Third, report the PCM ratio in the supplementary material. A ratio near 1 provides evidence that the chains have explored the full posterior support. A ratio below 0.5 warrants investigation.

5.4 Limitations

First, our synthetic targets are Gaussian mixtures with known component parameters. Real multimodal posteriors have mode shapes and separations that are unknown a priori, and the PCM's effectiveness may differ for non-Gaussian modes with heavy tails or skewness.

Second, the PCM relies on the minimum enclosing ellipsoid, which is sensitive to outliers and may overestimate volume in heavy-tailed distributions. Robust covariance estimators (e.g., minimum covariance determinant) could improve the PCM at the cost of computational complexity.

Third, the reference volume $V_{\text{prior}}$ requires a prior sample, which may be expensive to generate for complex models. In hierarchical models, sampling from the prior itself can be non-trivial.

Fourth, the $\text{PCM}_{\text{ratio}} < 0.5$ threshold was calibrated on Gaussian mixtures and may need adjustment for other target families. Adaptive thresholds based on dimension and sample size would be preferable but require further simulation study.

Fifth, the PCM does not identify which modes are missing or how many — it only flags that the posterior volume is smaller than expected. Follow-up analysis with tempering or nested sampling is needed to characterize the missing modes.

6. Conclusion

Standard MCMC convergence diagnostics produce an 18% false convergence rate on multimodal targets with three or more modes, rising to 29% for eight-mode targets in 50 dimensions. Split-$\hat{R}$ reduces the rate to 11% but remains inadequate for high-dimensional multimodal targets. The Posterior Contraction Monitor, which tracks the volume expansion of posterior samples, achieves a 4.8% false convergence rate with modest computational overhead. Combining the PCM with split-$\hat{R}$ brings the rate below 3%. These results argue against using any single convergence diagnostic and in favor of a multi-criterion approach that includes at least one volume-based or support-based check when multimodality is plausible.

An implementation of the PCM for Stan, PyMC, and NumPyro is available at https://github.com/spike-tyke/pcm.

References

[1] Gelman, A. and Rubin, D.B., 'Inference from Iterative Simulation Using Multiple Sequences,' Statistical Science, 1992.

[2] Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C., 'Rank-Normalization, Folding, and Localization: An Improved $\hat{R}$ for Assessing Convergence of MCMC,' Bayesian Analysis, 2021.

[3] Geweke, J., 'Evaluating the Accuracy of Sampling-Based Approaches to the Calculation of Posterior Moments,' in Bayesian Statistics 4, Oxford University Press, 1992.

[4] Brooks, S.P. and Gelman, A., 'General Methods for Monitoring Convergence of Iterative Simulations,' Journal of Computational and Graphical Statistics, 1998.

[5] Skilling, J., 'Nested Sampling for General Bayesian Computation,' Bayesian Analysis, 2006.

[6] Yao, Y., Vehtari, A., Simpson, D., and Gelman, A., 'Stacking of Predictive Distributions,' Bayesian Analysis, 2018.

[7] Margossian, C.C., Hoffman, M.D., Sountsov, P., Riou-Durand, L., Vehtari, A., and Blei, D.M., 'Nested $\hat{R}$: Assessing the Convergence of Markov Chain Monte Carlo When Running Many Short Chains,' Bayesian Analysis, 2024.

[8] Carpenter, B. et al., 'Stan: A Probabilistic Programming Language,' Journal of Statistical Software, 2017.

[9] Neal, R.M., 'MCMC Using Hamiltonian Dynamics,' in Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, 2011.

[10] Hoffman, M.D. and Gelman, A., 'The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo,' Journal of Machine Learning Research, 2014.

[11] Robert, C.P. and Casella, G., Monte Carlo Statistical Methods, 2nd ed., Springer, 2004.

[12] Betancourt, M., 'A Conceptual Introduction to Hamiltonian Monte Carlo,' arXiv:1701.02434, 2017.