When Sampling Fails

Matthew Reda Bayesian workflow

Every applied Bayesian eventually stares at a fit that refuses to converge: R-hat stuck at 1.4, a scatter of divergences, an effective sample size in the double digits. The instinctive responses are to run longer, add chains, and raise target_accept until the warnings stop. All three are usually wrong as a first move. The diagnostics are not interchangeable alarm bells: each one points at a different disease, and only one of the three is ever cured by more iterations. This post reads them properly, then walks the escalation ladder from cheapest fix to last resort. It ends with Sequential Monte Carlo, the exact sampler you call for a second opinion rather than for speed.

The Sampler Is an Instrument, and It Is Telling You Something

The organizing idea comes from what Gelman and colleagues (2020) call the folk theorem of statistical computing: when you have computational problems, often there is a problem with your model. A sampler is a measurement instrument pointed at a posterior distribution. When the instrument struggles, the overwhelmingly common cause is a posterior whose geometry is pathological: a funnel it cannot descend into, a ridge it can only crawl along, or two separated modes it cannot travel between. The instrument itself is rarely the problem. Modern implementations of the No-U-Turn Sampler (Hoffman & Gelman, 2014) are extraordinarily well engineered.

That reframing has a practical consequence. A failed fit looks like an obstacle between you and your results. Usually it is a result. It tells you that the model, as parameterized, asks for something the data cannot deliver, and the specific diagnostic that fired tells you which kind of something it asked for. The escalation ladder below is ordered accordingly: understand the signal first, spend seconds before minutes, and fix the model's identification before reaching for a different algorithm.

Three Numbers, Read Correctly

Modern convergence assessment rests on three quantities, standardized in their current form by Vehtari and colleagues (2021).

Definition: Split R-hat

Run \( M \) chains of \( N \) draws, split each in half (so slow drift within a chain shows up as disagreement), and compare the between-chain variance \( B \) to the within-chain variance \( W \):

$$ \hat{R} \;=\; \sqrt{\frac{\widehat{\mathrm{var}}^{+}}{W}}, \qquad \widehat{\mathrm{var}}^{+} \;=\; \frac{N-1}{N}\,W \;+\; \frac{1}{N}\,B . $$

If every chain explores the same distribution, \( B \approx W \) and \( \hat R \to 1 \). Vehtari and colleagues (2021) compute it on rank-normalized values, a form of the statistic that holds up under heavy tails and infinite variance, and they recommend acting only when \( \hat R \lt 1.01 \). That bar is deliberately stricter than the older 1.1 convention, which they show can pass fits whose estimates are still badly unstable.

Effective sample size answers a different question: of the \( MN \) draws you paid for, how many independent draws are they worth? Autocorrelation \( \rho_t \) at lag \( t \) discounts them,

$$ \mathrm{ESS} \;=\; \frac{MN}{\,1 + 2\sum_{t=1}^{\infty} \rho_t\,} , $$

and the practical floor is roughly 400 effective draws in the bulk (about 100 per chain at four chains) before R-hat itself is even reliably estimable. Check tail-ESS separately, because the quantiles that make up your credible intervals converge more slowly than the mean. Low ESS with good R-hat gets reported as failure to converge. What it shows is slow mixing. The chains agree about where the posterior is, and they are simply expensive to move around in.

Divergences are the odd one out: an event counter rather than a statistic, and specific to Hamiltonian samplers. HMC simulates a physical trajectory through the posterior. When the true curvature changes faster than the integrator's step size can track, the simulated energy explodes and the trajectory is flagged as divergent (Betancourt, 2017). What matters is where divergences happen. They cluster exactly in the regions the sampler is being repelled from: the neck of a funnel, the corner of a constraint. A divergent fit is therefore biased precisely in the part of parameter space you most needed explored, and no amount of additional sampling fixes that, because the additional samples avoid the same region.

Three Failure Syndromes

Put the three numbers together and the differential diagnosis falls out.

PresentationLikely diseaseFirst-line treatment
Divergences, often with high R-hatPathological curvature: funnels in hierarchical variances, cliffs from constraints or mis-scaled dataReparameterize (non-centered), rescale. target_accept only after
High R-hat, but each chain individually mixes cleanlyMultimodality: each chain found a different mode (label switching, sign reflections, aliased structure)Identification (anchors, constraints, ordering), then confirm with SMC
Low ESS across many parameters, R-hat roughly fineA ridge: two parameters trading off (weak identification), or strong posterior correlationPriors that break the ridge, more information, then and only then longer runs

The middle row is where the standard advice actively misleads. When two chains sit in two modes, each one looks healthy in isolation: good within-chain mixing, stable traces, no divergences. Only the between-chain comparison catches it, which is precisely why R-hat is computed across multiple independently-initialized chains in the first place. Doubling the draws does nothing. The chains simply converge ever more confidently to their separate answers. In mixture and latent-structure models the situation is generic rather than exotic. Permute the class labels, or flip the sign of a latent factor and its loadings, and the likelihood is exactly unchanged (Stephens, 2000). The posterior genuinely has \( k! \) or \( 2^k \) copies of every mode, and a gradient-based chain started near one copy will stay there.

Two four-hundred-iteration traces, two chains each. Left: a healthy fit. Both chains explore the same region, so between-chain variance matches within-chain variance and \( \hat R \approx 1.00 \). Right: textbook multimodality. Each chain mixes cleanly on its own but sits in a different mode, so \( \hat R \approx 1.9 \) even though no single trace looks pathological. Only the between-chain comparison catches it, which is exactly why \( \hat R \) is computed across independently-initialized chains.

💡 The wall-clock diagnostic

There is a fourth signal the diagnostics tables never list: time. A model whose warmup will not settle, or whose sampling pace is an order of magnitude slower than a comparable model, is exhibiting the same geometry problems before any threshold trips. Treat "this is inexplicably slow" as a divergence-class symptom rather than an infrastructure complaint. The folk theorem applies to runtime too.

Cheap Reconnaissance: MAP, Laplace, and Variational Fits

Before touching the sampler's dials, it pays to interrogate the posterior with tools that cost seconds. None of them produce calibrated uncertainty (in the framework they are all flagged approximate=True, and the reports refuse to present them as final inference), but each one answers a specific diagnostic question.

A MAP optimization asks the cheapest question of all: does a well-behaved mode even exist, and can gradient ascent find it? An optimizer that stalls or returns NaN almost always indicates a mechanical defect, such as data on wild scales, a gradient that is undefined at the initialization, or a constraint boundary. Any of those would have cost an hour of NUTS warmup to discover the slow way.

A Laplace approximation (Tierney & Kadane, 1986) goes one step further: expand the log posterior to second order around the MAP point \( \hat\theta \), giving

$$ p(\theta \mid y) \;\approx\; \mathcal{N}\!\left(\hat\theta,\; H^{-1}\right), \qquad H \;=\; -\nabla^2 \log p(\theta \mid y)\big|_{\hat\theta} . $$

The Gaussian is a caricature, but the attempt is informative. If \( H \) fails to be positive definite, some direction of parameter space is flat at the mode. That is a weak-identification flag delivered in seconds, naming the parameters involved, before any MCMC has run. A Laplace fit is also a steadier quick-check than bare MAP on high-dimensional models, because the curvature draws smooth over the optimizer's run-to-run jitter.

Two variational fits give a cheap sketch of the posterior's shape and location: ADVI (Kucukelbir et al., 2017) and the quasi-Newton Pathfinder (Zhang et al., 2022). An ELBO that oscillates without ever improving is another geometry smell. A variational posterior that lands somewhere absurd suggests the priors and the likelihood are fighting each other. In the framework all four run through the same interface:

# illustrative
results = model.fit(method="map")       # does a mode exist?
results = model.fit(method="laplace")   # is the curvature sane? (Hessian = identification probe)
results = model.fit(method="advi")      # cheap posterior sketch
results.approximate                     # True: never final inference

Fix Identification, Not the Sampler

The rung where failed fits actually get fixed is the model itself, and the fixes are boringly repeatable. The canonical example is the hierarchical funnel. In a centered parameterization, group effects \( \theta_g \sim \mathcal{N}(\mu, \tau) \) create a geometry where the plausible scale of every \( \theta_g \) depends on \( \tau \): as \( \tau \to 0 \) the posterior narrows into a funnel neck whose curvature defeats any fixed step size (Betancourt & Girolami, 2015). The non-centered form,

$$ \theta_g \;=\; \mu + \tau \,\eta_g, \qquad \eta_g \sim \mathcal{N}(0, 1), $$

is algebraically the same model with the dependence between location and scale removed from the prior, and it routinely converts a divergence-riddled fit into a clean one. Which form samples better depends on how strongly the data constrain each group (Papaspiliopoulos, Roberts & Sköld, 2007): well-measured groups favor centered, weakly-measured groups favor non-centered. Either way the cure is a reparameterization, and no sampler setting substitutes for it.

Neal's funnel: why centered fails and non-centered fixes it

Each dot is a draw of \( (x, v) \) where \( v \sim \mathcal{N}(0, \sigma_v) \) and \( x \sim \mathcal{N}(0, e^{v/2}) \), the geometry a hierarchical variance creates. Deepen the funnel and its neck (bottom) pinches to a width no fixed integrator step can resolve: in the centered parameterization the sampler is flung out and diverges there, biasing exactly the region it can't enter. Toggle to non-centered \( x = e^{v/2}\eta \) and you get the same posterior with regularized geometry, divergences gone.

3.00
Parameterization…
Divergences…
Neck width (x-scale at base)…
Sampler verdict…

Divergences cluster where the sampler is repelled, not at random, so more draws only avoid the same neck. The fix is the reparameterization, not target_accept.

Multimodality has its own family of identification fixes. Label switching in a mixture is removed by an ordering constraint or a post-hoc relabeling rule (Stephens, 2000). Sign reflection in a latent factor, the \( (F, \lambda) \leftrightarrow (-F, -\lambda) \) symmetry, is removed by anchoring one loading positive. And the anchor has to sit where the data actually pin that loading down, which is a stricter requirement than it first sounds. In our own structural mediation work, a factor anchored on a loading the data held near zero still escaped into its reflected mode, because a positivity constraint costs nothing at zero. That showed up as a split R-hat of 1.755 with individually clean chains. Re-anchoring on a loading the data held firmly nonzero dropped R-hat to 1.08. Ridges respond to the same logic. When two parameters can trade off (a carryover rate against an autoregressive coefficient, a scale against a loading), the fix is a prior or constraint that breaks the trade, or an experiment that informs one side of it.

⚠️ A different sampler does not remove modes

It is tempting to treat multimodality as a sampler deficiency: NUTS gets stuck, so find an algorithm that doesn't. But if the modes come from a symmetry of the model, they are exact copies. Sampling all of them perfectly yields a posterior mean of zero for every symmetric quantity, which is cancellation dressed up as inference. The sampler that reaches every mode (below: SMC) is the diagnosis. An identification constraint is the treatment. Only when modes are semantically real, genuinely competing explanations the model should entertain, is the multimodal posterior itself the answer.

When Sampler Settings Are the Answer

Settings have a legitimate but narrow role, and the order matters. Raising target_accept (0.9 → 0.95 → 0.99) shrinks the step size so the integrator can track tighter curvature. That is the right move for a few stubborn divergences on an otherwise-healthy fit, once reparameterization has done the heavy lifting. More warmup helps a sampler whose adaptation had not settled. Extra draws and chains are the honest cure for a pure ESS shortfall on a slow-mixing but correct geometry. A faster backend (the framework defaults to NumPyro's NUTS) buys wall-clock, not correctness.

The failure mode to avoid is using settings to silence rather than cure. On an unfixed funnel, target_accept=0.999 trades divergence warnings for a sampler that crawls, and the bias that remains is now invisible. If a setting change makes the warnings disappear, go back and re-run the diagnostics that were failing before, because those are the ones that confirm the fix. An absent warning confirms nothing on its own.

Sequential Monte Carlo: The Second Opinion

When the syndrome table says "multimodality" and the identification fixes are in dispute, you want an exact method that cannot get mode-locked. That is what tempered Sequential Monte Carlo is for (Del Moral, Doucet & Jasra, 2006; Chopin, 2002; the tempering idea traces to Neal's 2001 annealed importance sampling). Instead of one point moving through the posterior, SMC evolves a whole population of particles through a sequence of distributions that interpolate from the prior to the posterior,

$$ p_{\beta}(\theta) \;\propto\; p(\theta)\; p(y \mid \theta)^{\beta}, \qquad 0 = \beta_0 \lt \beta_1 \lt \cdots \lt \beta_T = 1 . $$

At \( \beta = 0 \) the particles blanket the prior, so every mode's basin is populated before the modes exist. As \( \beta \) rises, the likelihood switches on gradually. Importance weighting, resampling, and short MCMC moves let the population concentrate into whatever structure the posterior has, with mass allocated across modes in proportion to what they actually contain. A gradient chain must travel between modes through a probability desert. An SMC population never has to make the trip.

The framework runs it as fit(method="smc") and launches several independent SMC runs, computing R-hat across them. SMC is exact, so nothing about the result is flagged approximate. That cross-run R-hat is the second opinion. If independent tempered populations agree with each other but not with your NUTS chains, NUTS was stuck. When the SMC runs also disagree, the posterior's multimodality is severe enough that even population methods struggle, and the model needs identification work regardless of sampler.

# illustrative
results = model.fit(method="smc", draws=2000, chains=4)   # 4 independent runs x 2000 particles
results.approximate                                       # False: exact sampler
results.diagnostics["rhat_max"]                           # across runs: the multimodality verdict
results.diagnostics["log_marginal_likelihood"]            # model evidence (see below)

Honesty about the costs. SMC is not a speedup and should never be chosen as one. On a well-behaved unimodal posterior, NUTS is far more sample-efficient, because gradients are simply the best guidance system available. SMC's particle moves (an independent Metropolis-Hastings kernel in the PyMC implementation) degrade with dimension, so a model carrying many hundreds of free parameters (a large geo hierarchy, a latent-state structural model) will crawl or worse. Its niche is exactly the two things gradient MCMC cannot do: populate separated modes, and estimate the quantity in the next section.

Model Evidence and Bayes Factors

The tempering ladder yields a by-product that no other sampler in the workflow provides. Each increment \( \beta_t \to \beta_{t+1} \) estimates a ratio of normalizing constants from its importance weights, and the ratios telescope into the marginal likelihood, the normalizing constant of Bayes' theorem itself:

$$ Z \;=\; p(y) \;=\; \int p(y \mid \theta)\, p(\theta)\, d\theta, \qquad \frac{Z_1}{Z_0} \;=\; \prod_{t=0}^{T-1} \frac{Z_{\beta_{t+1}}}{Z_{\beta_t}} . $$

The ratio of two models' marginal likelihoods is the Bayes factor (Kass & Raftery, 1995), the factor by which the data shift the odds between them, with the integral's built-in penalty for unused flexibility. In a failed-fit investigation this is the tool for adjudicating between candidate fixes: two adstock families, two mediation structures, a one- versus two-mode specification. Marginal likelihoods are famously sensitive to priors, since a diffuse prior on an extra parameter is punished whether or not the data use it. Treat a Bayes factor from SMC as evidence alongside predictive comparison, never as a verdict. But when the question is structural ("does the second mode earn its place?"), the evidence integral is the right instrument, and SMC computes it as a side effect of sampling.

Where an MMM Breaks, Specifically

A media mix model assembles most of the classic pathologies in one place, which is why it makes a good closing case study.

Ridges are the default here. Adstock carryover and an autoregressive baseline both say "this week resembles last week," so their parameters trade off along a ridge the data barely constrain, which surfaces as low ESS with R-hat sitting near the line. A saturation half-point and a coefficient trade off similarly when spend never leaves a narrow band. The durable fixes are informative: an ROI-scale prior, a lift-test calibration that pins one end of the trade (see the calibration post), or spend variation designed into the media plan itself.

Hierarchy means funnels. Geo-level pooling puts the eight-schools geometry at the center of the model. Per-geo effect scales are exactly the \( \tau \) whose funnel neck generates divergences. Non-centered parameterization is the standard cure, and the framework's structural models choose parameterization by how densely each latent state is measured.

Latent structure means symmetry. Mediators, latent demand factors, and latent-class extensions bring reflection and relabeling modes, which is the high-R-hat-with-clean-chains syndrome from the middle row of the table. Sign anchors placed where the data pin a loading down are the treatment. SMC across independent runs is the confirmation.

The pre-fit checks deserve the last word, because they are where much of this post becomes unnecessary. Prior predictive simulation and simulation-based calibration (Talts et al., 2018) catch a large share of "sampling problems" before the first real fit, whether the trouble is priors implying absurd data or a computation that was never calibrated in the first place. A companion post covers both in depth. The framework's troubleshooting guide encodes the full ladder (diagnostics, approximate reconnaissance, identification, settings, SMC) as its standard operating procedure for a fit that will not converge.

Takeaways

  • The three diagnostics are not interchangeable: divergences mean broken curvature, high R-hat with clean chains means multimodality, low ESS means a ridge. Only the ESS shortfall is ever cured by running longer.
  • Divergences cluster where the sampler is repelled, so a divergent fit is biased exactly where it matters most, and more draws only avoid the same region.
  • Reconnaissance is nearly free: MAP asks whether a mode exists, a Laplace fit's failing Hessian names the weakly-identified directions, variational fits sketch the shape. None are final inference.
  • The durable fixes are identification: non-centered parameterizations for funnels, sign anchors placed where the data pin a loading down, priors or experiments that break ridges. Settings like target_accept are the last dial, not the first.
  • SMC is the second opinion, not a speedup. Independent tempered particle populations cannot get mode-locked, so R-hat across SMC runs separates "NUTS was stuck" from "the posterior really has modes." The tempering ladder also yields the marginal likelihood for Bayes-factor comparison of candidate fixes.
  • An MMM concentrates all three syndromes: adstock↔baseline ridges, geo-hierarchy funnels, and latent-structure symmetries. Pre-fit checks (prior predictive, SBC) catch a large share before the first real fit.

References