The Simulation Loop Does More Than SBC

Matthew Reda Bayesian workflow

Simulation-based calibration costs one full refit per simulated dataset. This framework’s run_mmm_sbc defaults to 64 of them, and on a national media mix model those refits run for hours. What comes back is a single verdict: the sampler is calibrated, or it is not, averaged over datasets drawn from the model’s own prior. Three other jobs run on that same loop. Two of them come after the fit: one conditions the calibration on the dataset you actually have, and one turns the ranks into a correction you apply to a cheap posterior instead of a verdict you can only obey. The third runs before any data exist at all, and asks whether the model could ever answer the budget question.

The Ensemble You Are Grading

Start with what a prior-SBC ensemble is made of. In the framework’s ROI parameterization, which the agent’s spec-built models select, each channel carries roi_<ch> ~ LogNormal(0, 1) (config/model.py:191, 194-195), a defensible weakly-informative choice on the decision scale. Draw 20,000 worlds from it: eight channels, 104 weeks, media spend at 35% of a revenue index whose baseline sits at 100, a normal baseline prior centered there with a standard deviation of 30, and half-normal observation noise. The implied R² runs from 0.18 at the fifth percentile to 1.00 at the ninety-fifth, with a median of 0.80, so these are recognizably weekly sales series.

Then look at the businesses they describe. In 69% of the worlds at least one channel returns more than three dollars per dollar. Better than five dollars per dollar in 36% of them. The median world’s best channel sits at 4.0x, and the ninety-ninth percentile at 20.6x. Those are not media plans anyone has run.

The other half of the ensemble is worse, because it is unanswerable rather than implausible. In 74% of the worlds some channel’s weekly contribution has a standard deviation smaller than a tenth of the observation noise, which is to say the channel is invisible in its own simulated data. Ask for both conditions at once, no channel above 3x and every channel at least a tenth of the noise floor, and 7.4% of the ensemble survives. Of the 64 refits run_mmm_sbc pays for, roughly five land on a world an analyst would recognize as their client’s.

Where this lives in the code

diagnostics/sbc.py implements Talts-style prior SBC and nothing else: run_mmm_sbc takes a built, unfitted BayesianMMM, calls sample_prior_predictive once for all simulations, then swaps the observed series with pm.observe on the same graph and refits (defaults n_sims=64, L=100). Everything below is about what else that loop could be doing.

SBC is doing exactly what it was built to do, and the priors are fine. Talts and colleagues designed the method for people writing samplers, and for them averaging over the prior is the right thing: an algorithm that only works on the plausible corner of parameter space is an algorithm with a hole in it. Säilynoja, Schmitt, Bürkner and Vehtari (2025) draw the line where it belongs. Prior SBC is a tool for algorithm developers, posterior SBC is a tool for modelers, and there are orders of magnitude more modelers than there are people writing NUTS implementations. The rank statistic and the histogram shapes belong to the companion post. Take them as given here.

The implausible worlds are also, and not coincidentally, the slow ones. On a Lotka-Volterra model Säilynoja and colleagues report 250 prior-SBC iterations taking 6.5 hours against 500 posterior-SBC iterations taking 2.5 hours. The weakly-informative prior puts real mass on parameter values whose posteriors are pathological to sample, so the sampler grinds. Conditioning on data removes that mass, and their runtime per iteration came out more than five times faster.

Conditioning on the Data You Have

The reformulation is short enough to state in four lines, and the reason it works is sequential updating: yesterday’s posterior is today’s prior, so the uniformity argument that underwrites SBC applies with \( \pi(\theta \mid y_{\text{obs}}) \) standing in for the prior.

Definition: Posterior SBC

Given an observed dataset \( y_{\text{obs}} \) and a fitted posterior \( \pi(\theta \mid y_{\text{obs}}) \), repeat:

$$ \theta' \sim \pi(\theta \mid y_{\text{obs}}), \qquad y_i \sim \pi(y \mid \theta'), \qquad \theta \mid (y_i, y_{\text{obs}}) \;\text{refit}, \qquad r_i = \operatorname{rank}\big(f(\theta');\, f(\theta_1), \dots, f(\theta_L)\big). $$

The rank of the drawn \( \theta' \) among draws from the augmented posterior is uniform when the computation is correct, for the same reason the prior-SBC rank is. Conditioning changes what gets graded. The sampler now sees datasets that resemble the one you have, because \( \theta' \) came from a posterior that has already seen it.

What that buys is sensitivity. Take two regions of the prior whose posteriors are biased in opposite directions. Prior SBC averages over both, the errors cancel in the rank histogram, and the histogram stays flat however far you push the iteration count, because more iterations estimate the average more precisely, and the average they are estimating is zero. Posterior SBC visits one region, so nothing cancels. A rank histogram is an average, and averages hide sign changes.

Their hierarchical case study shows the same blindness arriving by dilution rather than cancellation, and it is the one to show a geo modeler. After 500 iterations of prior SBC, Säilynoja and colleagues report no noticeable calibration issue with either the centered or the non-centered parameterization, every PIT-ECDF difference line inside the simultaneous band, which would tell the modeler both implementations are fine. Posterior SBC conditional on a dataset with a weak population prior and a strong likelihood then finds the centered version calibrated and the non-centered one failing on the population-level parameters. Conditional on the other dataset, strong prior and weak likelihood, the verdict reverses.

Which parameterization is right is therefore a fact about the data, settled after the fit rather than before it. A geo MMM with per-geo channel effects is this model, and the information per geo is a property of the panel you were handed. The sampling-failure post treats non-centering as a repair you reach for when the diagnostics complain. Posterior SBC lets you pick the parameterization your panel supports before they have to, and it is the only check in the toolkit that can, because it is the only one that conditions on the panel.

# illustrative
# Posterior SBC on top of an already-fitted model.
post = fit(y_obs)                                   # the fit you were going to ship
ranks = []
for i in range(n_sims):
    theta_prime = post.draw()                       # NOT a prior draw
    y_i         = simulate_likelihood(theta_prime)  # a dataset like the one you have
    aug         = fit(concat(y_obs, y_i))           # refit on the augmented data
    draws       = aug.draws(quantity=joint_loglik)  # test quantity, not just params
    ranks.append(sum(d < joint_loglik(theta_prime) for d in draws))

Note the test quantity in that last line. The 2025 paper uses the joint log-likelihood throughout as its default, which is the recommendation Modrák and colleagues arrived at from the other direction, and which the companion post covers. Those two results compose cleanly. Rank a data-dependent quantity, and rank it on data that resembles yours.

The framework already has the single-point version

run_recovery_coverage(model, truth="posterior_mean") in diagnostics/coverage.py fixes every free parameter at the fitted posterior mean with pm.do, simulates datasets there, refits each, and measures interval coverage. That is posterior SBC collapsed to one point of the posterior instead of draws from it, which makes it cheaper and blind to anything that varies across the posterior. Moving from a point to draws is the change, and the loop is otherwise the same code.

From Check to Correction

Every diagnostic above returns a verdict. Cai, Greengard, Goodrich and Gelman (2026) turn the same ranks into a repair for the posterior.

The move is a single scalar. Take an approximate posterior, rescale its draws around their own mean by a factor \( k \), and choose \( k \) so that the SBC quantiles come out nominal. Two selection rules work. Grid search directly on the coverage objective is the obvious one. The cheaper one needs no search at all: compute the SBC z-scores \( (\theta_{\text{true}} - \mu_{\text{post}}) / \sigma_{\text{post}} \) across replications and take their standard deviation, which is the factor by which the intervals are too narrow. On eight schools, with Stan’s ADVI on a centered parameterization and HMC on the non-centered version as ground truth, the two rules agree: grid search lands between 2.41 and 2.56, the z-score rule at 2.48. Under that rule, adjusted coverage comes out at 0.958, 0.909, 0.788 and 0.505 against nominal levels of 0.95, 0.90, 0.80 and 0.50.

What makes 2.4 interesting is that ADVI’s posteriors on that problem are about half as wide as HMC’s in their Figure 3a, and half as wide should want a factor of 2. The extra comes from location. ADVI’s posterior mean is shifted, and the shift is different in every replication, so a correction that can only change the width has to over-widen to swallow the wandering center. I wanted to see that mechanism on its own, so I planted both errors deliberately in a conjugate normal-normal model and ran the framework’s own run_sbc over 20,000 replications: an approximate posterior half as wide as the exact one, with a mean displaced each time by a normal draw of half the exact posterior’s standard deviation.

Nominal levelEmpirical coverage, broken posteriorAfter widening by 2.31
50%24.0% [23.4, 24.6]50.9%
80%44.0% [43.3, 44.7]80.0%
90%54.7% [54.0, 55.4]89.6%
95%63.9% [63.2, 64.5]94.2%

Brackets are Jeffreys binomial intervals from coverage_from_ranks. The 90% interval covering 54.7% is the exact symptom technical-docs/coverage-diagnostics.md is written to diagnose, and the diagnosis it gives, an approximate fit, is correct here by construction.

Three ways of picking \( k \) agree with each other and disagree with 2. Grid search on the coverage objective gives 2.31, the z-score rule gives 2.22, and the closed form for this planted setup, \( \sqrt{1 + 0.5^2} / 0.5 \), gives 2.236. Now rerun the identical loop with the location error switched off and the width error untouched. The z-score rule returns 1.99. Between that and 2.22 sits the price of the wandering center, stated in the units of the correction you have to apply. That is a more honest account of what a variational fit costs than “uncertainty is not calibrated,” which is what the framework’s own troubleshooting table currently says.

Deep diveWhy mean-field intervals are narrow by theorem

Margossian, Pillaud-Vivien and Saul (2025) prove an impossibility result for factorized approximations. If the target posterior does not factorize, no factorized \( q \) can match more than one of three uncertainty measures at a time: the marginal variances, the marginal precisions, or the generalized variance. Minimizing \( \mathrm{KL}(p \,\|\, q) \) matches the marginal variances. Minimizing \( \mathrm{KL}(q \,\|\, p) \), which is the standard variational objective and the one ADVI optimizes, matches the marginal precisions instead. A marginal precision is a conditional variance in disguise, so matching it and reporting it as a marginal is how the intervals come out too tight.

An MMM posterior is the textbook non-factorizing case. Channel spend series move together, adstock decay trades off against saturation, saturation trades off against the coefficient, and the joint sits on a ridge rather than a peak. The theorem therefore says the narrow ADVI credible intervals you get from an MMM are structural, and no amount of tuning the step size or the iteration count will widen them. The best response is to stop using mean-field variational inference wherever you intend to report uncertainty. Failing that, use full rank and pay for it. Failing that, use mean field and apply the correction above. The factor absorbs a location error as well as a width error, and only the width part is honestly a scale.

⚠️ The obvious extension of this is invalid

If SBC ranks averaged over the prior can calibrate a correction, the natural next thought is to average over the posterior instead and get a correction tuned to your own data. Cai and colleagues prove it fails. In a plain normal-normal model, swapping the prior for the posterior in the averaging step produces z-scores that are not standard normal for any value of \( y \), and recalibrating against them yields intervals that come out both shifted and overconfident. Worse, that same model quietly pools only half as far toward the prior as the exact posterior, its mean moving from \( y/2 \) to \( 3y/4 \). The recalibration has edited the inference rather than the interval drawn around it. Hierarchical models are the one case they think could survive. Posterior SBC is the move that validly conditions on your data, because it changes what gets simulated. This changes what the ranks mean.

The Question Before the Data

Betancourt’s principled workflow opens by asking four questions of a model, and this blog has answered two of them. Question One asks whether the model is consistent with domain expertise, and prior predictive checks answer it. Question Two asks whether the computational tools will be sufficient to fit the posterior accurately, and SBC answers it. Question Four asks whether the model is rich enough to capture the true data generating process, and posterior retrodictive checks answer that one later, on real data. Question Three has no coverage anywhere in this series: will our inferences provide enough information to answer our questions?

Betancourt answers it in Step Eleven of the workflow, inferential calibration, run over the simulated ensemble that Step Six already built for algorithmic calibration in Step Ten. The loop is the same loop. Simulate from the prior, fit each replicate, and instead of ranking the truth, measure how wide the posterior came out and whether that width is small enough to move a decision.

For a media mix model the decision is a budget shift, and a budget shift has a threshold. Distinguishing a channel that returns 1.2x from one that returns 3x is 1.8 ROI points. If the posterior interval on that channel’s ROI is wider than 1.8 no matter what the data say, the model has failed before anyone has collected anything, and it has failed in a way that no amount of sampler correctness will repair.

This one is computable in closed form for a linear-Gaussian model, which is the useful special case, because with a Gaussian prior and known noise the posterior covariance does not depend on \( y \) at all. The design: 104 weeks, eight channels at spend shares from 32% down to 3%, geometric adstock and a Hill curve, columns rescaled so the coefficient is the channel’s ROI in revenue per media dollar, media at 12% of revenue, spend partly chasing a shared demand driver, an intercept, a linear trend, two Fourier seasonality pairs, and noise set to R² = 0.85. Prior on each ROI: normal, centered at 1.5, standard deviation 1.0, whose central 80% interval is 2.56 points wide.

DesignMedian pairwise correlation among media columnsChannels whose 80% ROI interval clears the 1.8-point gap
Base: 2 years, media at 12% of revenue0.901 of 8 (TV alone)
Double the history to 4 years0.891 of 8
Raise media to 25% of revenue0.901 of 8
Halve the noise (R² = 0.95)0.902 of 8
Double the flighting variation0.662 of 8
Decouple spend from demand entirely0.142 of 8
Four of those changes at once0.677 of 8

Read the middle rows first. Every single intervention buys at most one additional channel, and doubling the length of the history buys nothing at all, which is the opposite of what a client asks for when the intervals come back too wide. The eighth channel, at 3% of spend, has a posterior standard deviation of 1.00 in the base design against a prior standard deviation of 1.00. Data move it by nothing measurable, so its interval is the prior’s interval and every ROI number reported for it restates an assumption. The carryover-and-shape post makes the same point about the transform parameters, from the identification side.

Now read the last row. Four years of history, media at a quarter of revenue, twice the flighting variation, and a cleaner series together resolve seven channels of the eight. The effects compound. Ask “how do we get tighter ROI intervals” and the honest answer is never one thing. That is a conversation worth having before the modeling contract is signed rather than after the readout, and Question Three is where it belongs.

The framework computes the adjacent quantity already. planning/ holds the EIG and EVOI machinery that prices what an experiment would teach you. Pointing the same arithmetic at the observational design, and running it over the ensemble the SBC loop already generates, is the missing piece.

Worlds the Model Cannot Represent

Everything above simulates from the model being graded. That is the source of its power and the boundary of its reach, and it is why a clean SBC run certifies machinery rather than judgment. Gelman, Vehtari, McElreath and colleagues give the complementary stage its own place in the workflow, and in the 2020 preprint it is §4.3, separate from fake-data simulation in §4.1 and SBC in §4.2. That stage is about design and bias rather than code correctness.

Their worked example is 500 students with a midterm and a final. Plant a treatment effect of exactly 10, then break the design: let assignment depend on the pre-test through \( \Pr(z=1) = \mathrm{logit}^{-1}((x-50)/10) \), and the raw comparison of treated against control returns −13.8 with a standard error of 1.5. Linear adjustment for the pre-test rescues it to 9.7 with a standard error of 0.8. Then make the true \( E(y \mid x, z) \) nonlinear and change nothing else: the same linear adjustment under the same unbalanced design degrades to 7.3 with a standard error of 0.9, while under a balanced design it still recovers 10.5 with a standard error of 0.8. Same estimator, same planted truth, 27% low because the assignment mechanism and the functional form conspired.

The analogue is exact and uncomfortable. Media spend responds to expected demand, which is the unbalanced assignment mechanism, and adstock composed with saturation is a functional form nobody gets exactly right. The simulation returns a number rather than a verdict: how wrong does this estimator get when this specific assumption is wrong by this much. Gelman and colleagues frame it as the interesting question in the case where the check is guaranteed to fail by construction, a \( t_4 \) truth fit with a normal likelihood, and their phrasing is the whole argument for the stage: how bad will these inferences be?

The violation worlds ship with the framework

synth/dgp.py and synth/dgp_geo.py generate scenarios labeled by which assumption they break and whether the planted truth is representable by the model at all. The holdout post uses them as an acceptance gate, a pass-or-fail check against an external answer key before a model ships. Used the §4.3 way they answer a different question: run the world your model provably cannot represent, and record the size of the error rather than its existence. A bias you have measured is a bias you can carry into the recommendation as a stated haircut. A bias you have only detected is a bias you argue about.

Where the Budget Goes

Given one overnight window of compute and a fitted model, spend it on posterior SBC with the joint log-likelihood as the test quantity. It is cheaper per iteration than the prior version, and it grades the fit you are about to ship. On a geo model it also tells you which parameterization your panel supports. Where the fit is approximate because the full one will not finish, run the loop anyway and keep the z-scores. Their standard deviation is the widening factor, and it costs one line. The alternative is publishing a 90% interval that covers half the time. Read that factor as diagnosis as well as repair, because the excess over the ratio of the widths, the 0.23 that separated 2.22 from 1.99 here, is location error, and location error is not something a wider interval makes correct.

Before the data exist, ask Question Three. The linear-Gaussian version is closed form and runs in seconds, and it names the channels the proposed design can never resolve. That answer is worth more than any diagnostic available afterwards, because afterwards the only remaining move is to widen the interval and apologize.

Takeaways

  • Under the framework’s ROI-mode LogNormal(0, 1) prior, 69% of prior-SBC worlds contain a channel returning more than 3x and 74% contain a channel invisible in its own data. About five of the default 64 refits land on a business anyone runs.
  • Posterior SBC draws \( \theta' \) from the fitted posterior, simulates, refits on the augmented data, and ranks. It catches biases that cancel in a prior-SBC histogram, and on Lotka-Volterra it ran 500 iterations in 2.5 hours against 6.5 hours for 250 prior-SBC iterations.
  • Whether a hierarchical model wants the centered or the non-centered parameterization is a function of the observed data, so a geo MMM cannot settle it before the fit. Posterior SBC is the check that conditions on the panel.
  • SBC ranks yield a scalar interval-widening factor. In a planted normal-normal test the broken 90% interval covered 54.7% and widening by 2.31 restored it to 89.6%; the factor exceeds the width ratio of 2 because the approximate posterior’s center wanders, and with the wander switched off the same rule returns 1.99.
  • Betancourt’s Question Three is answerable before data collection. In a realistic eight-channel design, only one channel’s 80% ROI interval is narrow enough to separate 1.2x from 3x, doubling the history changes nothing, and four changes applied together resolve seven.
  • The stage that prices a bias simulates from a world the model cannot represent. Gelman and colleagues’ own example moves an adjustment estimator from unbiased to 27% low by changing only the assignment mechanism and the functional form.

References