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 — run longer, add chains, raise target_accept until the warnings stop — are usually the wrong first moves, because the three 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 the diagnostics properly, then walks the escalation ladder from cheapest fix to last resort — ending with Sequential Monte Carlo, the exact sampler you call not for speed but for a second opinion.

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 not a bug in the instrument — modern implementations of the No-U-Turn Sampler (Hoffman & Gelman, 2014) are extraordinarily well engineered — but 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.

That reframing has a practical consequence. A failed fit is not an obstacle between you and your results; it usually is a result. It is telling 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. 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 (robust to heavy tails and infinite variance) and recommend acting only when \( \hat R \lt 1.01 \) — a deliberately stricter bar 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 — plus a separate tail-ESS check, because the quantiles that make up your credible intervals converge more slowly than the mean. Low ESS with good R-hat is not "failure to converge"; it is slow mixing. The chains agree about where the posterior is; they are just expensive to move around it.

Divergences are the odd one out — not a statistic but an event counter, 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). The crucial property 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 structureIdentification: 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 deserves the emphasis, because it is the one 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 converge ever more confidently to their separate answers. In mixture and latent-structure models this is generic, not 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, not 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 — data on wild scales, a gradient that is undefined at the initialization, a constraint boundary — that would have wasted an hour of NUTS warmup to discover.

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 of the posterior, but the attempt is informative: if \( H \) fails to be positive definite, some direction of parameter space is flat at the mode — a weak-identification flag delivered in seconds, naming the parameters involved, before any MCMC has run. This is also why a Laplace fit is a more stable quick-check than bare MAP on high-dimensional models: the curvature draws smooth over the optimizer's run-to-run jitter.

Variational fits — ADVI (Kucukelbir et al., 2017) and the quasi-Newton Pathfinder (Zhang et al., 2022) — give a cheap sketch of the posterior's shape and location. An ELBO that oscillates without improving is another geometry smell; a variational posterior that lands somewhere absurd suggests the priors and likelihood are fighting. 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. The point is that the cure was a reparameterization, not a sampler setting.

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 \) — same posterior, 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 — 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 be chosen where the data have leverage: in our own structural mediation work, a factor anchored on a loading the data held near zero still escaped into its reflected mode (the positivity constraint costs nothing at zero), showing 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 not inference, it is cancellation. The sampler that reaches every mode (below: SMC) is the diagnosis; the 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 — the right move for a few stubborn divergences on an otherwise-healthy fit, after reparameterization has done the heavy lifting. More warmup helps a sampler whose adaptation had not settled; more 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: target_accept=0.999 on an unfixed funnel trades divergence warnings for a sampler that crawls, and the bias that remains is now invisible. If a setting change makes warnings disappear, the diagnostics that should confirm the fix are the ones that were failing before — not merely the absence of the warning.

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 — 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") — an exact method, so nothing about the result is flagged approximate — and launches several independent SMC runs, computing R-hat across them. 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; if 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 — 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 with 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 — a diffuse prior on an extra parameter is punished whether or not the data use it — so treat Bayes factors from SMC as one line of evidence alongside predictive comparison, not 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, not the exception. 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 — low ESS, R-hat 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 — the high-R-hat-with-clean-chains syndrome. Sign anchors placed where the data have leverage are the treatment; SMC across independent runs is the confirmation.

The pre-fit checks — prior predictive simulation and simulation-based calibration (Talts et al., 2018) — deserve the last word, because they are where much of this post becomes unnecessary. A large share of "sampling problems" are priors implying absurd data or a computation that was never calibrated in the first place, and both are detectable before the first real fit; the companion post covers them 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 — more draws 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 have leverage, 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" — and the tempering ladder 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