Validating calibration: simulation-based calibration

What the instrument reads

diagnostics() scores one fit. It asks whether the machinery that produced that posterior behaved: whether the chains mixed, whether the outer grid covered the hyperparameter posterior, whether the inner Gaussian was a reasonable shape for the latent conditional. Those are verdicts on the procedure at one data set.

Calibration is a different question, and no single fit can answer it. An inference algorithm is calibrated when its posteriors are self-consistent across the whole generative model: simulate a truth, simulate data at that truth, fit, and the truth should land at a uniformly distributed position inside the posterior it produced. Repeat that many times and those positions have a distribution you can test.

The usual way to check this is coverage: count how often a 95% interval contains the truth. That reads one point of the marginal CDF, or two if you also run a 50% interval. It cannot say whether a posterior is biased, over-dispersed, under-dispersed or asymmetric, and it is weak enough that two genuinely different reads of the same fit can score identically on it. sbc() reads the whole CDF.

The measurement is the probability integral transform. For a scored quantity whose true value is theta_0 and whose reported posterior CDF is F, take u = F(theta_0). Under exact inference u is exactly Uniform(0, 1), so the entire empirical CDF of the u values across simulations is the test statistic (Talts et al. 2018).

A model to validate

A balanced Gaussian random-intercept model: six regions of four observations each, an intercept and one covariate, a random intercept per region whose standard deviation sigma is drawn from a seven-point grid.

The design is small on purpose. The sigma posterior then spreads across four or five of the seven grid cells, which is the regime where the outer-grid mixture posterior for a fixed effect is farthest from the Gaussian matching its first two moments.

GRID <- exp(seq(log(0.2), log(1.5), length.out = 7))
PHI  <- 0.7
BETA <- c(-0.2, 0.7)

simulate_one <- function(seed) {
  set.seed(seed)
  sigma  <- GRID[sample.int(length(GRID), 1L)]
  region <- rep(seq_len(6L), each = 4L)
  X      <- cbind(1, rnorm(24L))
  u      <- rnorm(6L, 0, sigma)
  list(y      = as.numeric(X %*% BETA) + u[region] + rnorm(24L, 0, PHI),
       X      = X,
       region = region,
       theta  = c(beta1 = BETA[1], beta2 = BETA[2], sigma = sigma))
}

fit_one <- function(d) {
  tulpa_nested_laplace(
    y = d$y, n_trials = rep(1L, length(d$y)), X = d$X,
    prior = list(list(type = "iid", obs_idx = d$region,
                      n_units = max(d$region), sigma_grid = GRID)),
    family = "gaussian", phi = PHI,
    control = list(n_threads = 1L, keep_grid_hessians = TRUE,
                   auto_recenter = FALSE, progress = FALSE))
}

simulate_one() has to be a pure function of its seed, and it returns theta, the named vector of true values the PIT is taken against. auto_recenter = FALSE is not a speed knob: the argument that the sigma PIT is uniform needs the fitted grid to be the same seven points the truth was drawn from, and a recentred grid breaks it.

Reporting a fit in a shape sbc() can read

A fitter callback returns a named list of arms, each a named list over scored quantities, each entry one of the predictive shapes in ?sbc_predictive. Everything downstream dispatches on the shape’s kind tag, so a backend reporting something new is one shape rather than a parallel scorer.

Two of them cover a nested-Laplace fit. What the outer grid defines for a fixed effect is a Gaussian mixture over the cells, and tulpa_posterior_draws() realizes it: each draw picks a cell by its weight, then samples that cell’s inner Gaussian. coef() and vcov() report that same mixture’s first two moments, so building an sbc_normal() from them is the collapsed read of the identical posterior. For the hyperparameter the grid defines a distribution on a finite support, which is sbc_discrete().

arms <- function(d) {
  fit <- fit_one(d)
  m   <- coef(fit)
  se  <- sqrt(diag(vcov(fit)))
  D   <- tulpa_posterior_draws(fit, n = 2000)
  w   <- fit$weights / sum(fit$weights)

  list(
    mixture = list(
      beta1 = sbc_draws(D[, 1]),
      beta2 = sbc_draws(D[, 2]),
      sigma = sbc_discrete(as.numeric(fit$theta_grid), w)),
    collapsed = list(
      beta1 = sbc_normal(m[1], se[1]),
      beta2 = sbc_normal(m[2], se[2])),
    narrow = list(
      beta1 = sbc_normal(m[1], se[1] / 1.25),
      beta2 = sbc_normal(m[2], se[2] / 1.25)))
}

Every arm reads off one solve per seed, so an arm-to-arm difference carries no fit-to-fit noise. narrow is a deliberately mis-scaled posterior, the same moments with the standard deviation divided by 1.25. A calibration harness that cannot fail is worth nothing, so the run below has to put it outside the band.

The prior-predictive experiment

res <- sbc("prior_predictive",
           simulator  = simulate_one,
           fitter     = arms,
           n_sim      = 200L,
           flat_prior = c("beta1", "beta2"))
res

flat_prior is a guard doing its job. Ordinary SBC draws the truth from the prior, so an improper prior has nothing to draw from, and the nested-Laplace door puts no prior on the fixed effects. sbc() probes the simulator, sees that beta1 and beta2 do not move across simulations, and refuses to score them unless the caller names them. Naming them asserts that their flat prior leaves the PIT uniform by a structural argument; the assertion is checked in both directions and travels on the result in res$premises.

Three reads of the same PIT sample

ks and p_unif read the raw PIT ECDF. inside_folded and p_folded read 2 |u - 1/2|, which is also uniform under correct inference. The fold is where a symmetric error shows: a posterior that is too narrow pushes PIT mass towards both ends at once, which cancels in the raw ECDF and accumulates in the folded one. crps is the continuous ranked probability score, closed form for a Gaussian mixture, so the nested tier’s own posterior is scored with no Monte Carlo.

inside is the verdict against a simultaneous band. A pointwise binomial band is not one: at n = 100, holding each order statistic inside its own 95% interval holds all of them together only 44.71% of the time. The band here is calibrated by bisection against the exact crossing probability of the uniform order statistics, so an ECDF excursion anywhere along the curve is a 0.05-level event.

Discrete quantities are randomized within their atom, u = F(theta^-) + V P(theta), so sigma on its seven-point grid and a continuous coefficient share one uniform reference and one band. Reading a rank against a continuous uniform is the classic silent SBC bug, and it is not something you can opt into here.

The proper score, paired

summary(baseline = ) pairs the CRPS seed by seed against one arm.

summary(res, baseline = "mixture")

delta with its t is the proper-score verdict; a negative delta is the arm scoring better. p_sign is more powerful at detecting that two arms differ at all, but the sign test is not a proper score, so it cannot rank them.

Expect the mixture and collapsed arms to score close together. They carry the same first two moments, and the CRPS integrates the whole squared CDF difference, which those moments dominate. Where they separate is on the ECDF reads above, which is why both instruments are worth running.

The CRPS is a proper posterior score here only because the truth is drawn afresh each simulation. Hold the truth fixed across seeds and the CRPS-optimal forecast is a point mass at it, so a sharper wrong posterior wins. That is enforced rather than documented: a fixed-truth sweep is not offered as an experiment, and summary(baseline = ) refuses to rank one.

The picture

plot(res, arm = c("mixture", "narrow"), quantity = "beta2")

plot(res, arm = c("mixture", "narrow"), quantity = "beta2", folded = TRUE)

The plot draws the ECDF difference from uniform, so a calibrated read is a flat line at zero inside the band. The under-dispersed arm bows away from it, and the bow is larger folded than raw.

Calibration conditional on an observed data set

The experiment above averages over the prior. A user fitting their own data asks something narrower: is the inference reliable in the posterior geometry this data set produces. The prior average can miss a defect confined to a small region of parameter space, and it can flag one the observed data rules out.

experiment = "posterior" answers the narrow question (Sailynoja et al. 2026, Algorithm 2). Draw theta' from pi(theta | y_obs), simulate a replicate at theta', and take the PIT under the augmented posterior pi(theta | y, y_obs). That is ordinary SBC with pi(theta | y_obs) in the role of the prior, so the same band, the same folded read and the same proper score all carry over, and it needs no proper prior at all.

It takes six callbacks.

d_obs <- simulate_one(99L)

model <- list(
  data_obs = d_obs,

  fit = function(data) fit_one(data),

  draw_theta = function(fit, seed) {
    set.seed(seed)
    b <- tulpa_posterior_draws(fit, n = 1L)
    k <- attr(b, "cells")[1]
    c(beta1 = unname(b[1, 1]), beta2 = unname(b[1, 2]),
      sigma = as.numeric(fit$theta_grid)[k])
  },

  simulate = function(theta, seed) {
    set.seed(seed)
    region <- rep(seq_len(6L), each = 4L)
    X      <- cbind(1, rnorm(24L))
    u      <- rnorm(6L, 0, theta[["sigma"]])
    list(y = as.numeric(X %*% theta[c("beta1", "beta2")]) +
           u[region] + rnorm(24L, 0, PHI),
         X = X, region = region)
  },

  pool = function(obs, rep) list(
    y      = c(obs$y, rep$y),
    X      = rbind(obs$X, rep$X),
    region = as.integer(c(obs$region, rep$region + max(obs$region)))),

  arms = function(fit, data) {
    m  <- coef(fit)
    se <- sqrt(diag(vcov(fit)))
    D  <- tulpa_posterior_draws(fit, n = 2000)
    list(
      mixture = list(
        beta1 = sbc_draws(D[, 1]),
        beta2 = sbc_draws(D[, 2]),
        sigma = sbc_discrete(as.numeric(fit$theta_grid),
                             fit$weights / sum(fit$weights))),
      narrow = list(
        beta1 = sbc_normal(m[1], se[1] / 1.25),
        beta2 = sbc_normal(m[2], se[2] / 1.25)))
  },

  group_ids = function(data) data$region)

pres <- sbc("posterior", model = model, n_sim = 200L)
pres

The driver hands draw_theta and simulate different seeds, so set.seed(seed) at the top of each is the correct fixture. Giving both the same seed makes the replicate’s noise a function of the truth, which is not p(y | theta'), and it shows up as a non-uniform PIT with nothing wrong in the inference under test.

Two premises, and what checks them

Each premise silently turns the construction into something that is not SBC, so each has a guard, and each guard’s conclusion travels on the result.

str(pres$premises)

The augmented posterior conditions on both data sets. Fitting the replicate alone is ordinary SBC under a hand-made prior. sbc() refuses a pool() returning no more than the replicate, or no more than the observed data.

The replicate is conditionally independent of y_obs given theta. The nested tier integrates the random effects out, so theta carries no per-group value, and a replicate drawn on the same regions couples the two data sets through the group effects theta does not describe. simulate() above draws six fresh regions and pool() offsets their labels past the observed ones. Supply group_ids and the observable half of that is verified, the labels being disjoint. Omit it and the result records the premise as unverified rather than assumed. The other half, that the replicate’s group effects came from the prior rather than conditionally on y_obs, is not visible from outside the callback and is not claimed.

A third requirement is a property of draw_theta rather than a guard. theta' has to be a joint draw from pi(theta | y_obs). Sampling beta from tulpa_posterior_draws() and sigma independently from the grid weights gives the right two marginals and the wrong joint, and the sigma read leaves the band when you do it. Reading attr(b, "cells") takes the coefficient and the hyperparameter from the same cell, which is what makes the draw joint.

Where this sits beside the single-fit band

diagnostics() on a nested fit reports the outer Pareto-k-hat, the inner skewness estimate gamma_3 and the inner importance k-hat, combined into one reliability band. That band runs on the fit you already have; sbc() costs a few hundred fits. So the band is what you read routinely, and this is what you reach for when the answer has to hold up.

The cheap one is not a compressed version of the expensive one, and the two disagree in both directions. Measured over fifteen configurations at 1000 simulations (dev_notes/issue339/): a Poisson configuration carrying an outer k-hat of 0.196 with both inner scores in the good band, the cleanest verdict the band can give, fails calibration at p = 2.3e-13 on its intercept. A binomial configuration carrying an outer k-hat of 1.413, well past the 0.7 escalation threshold, passes at p = 0.17.

So read the shipped band as a screen. A rejection from sbc() is the strong statement: on that same measurement the false-positive rate is 0.0117 against a nominal 0.05, and power reaches 80% at roughly a 10% dispersion error or a 0.14-SD location bias, which makes a pass the weaker one.

diagnostics() reads both at once when you hand it the calibration result.

fit <- fit_one(d_obs)
diagnostics(fit, sbc = res)

See also

• ?sbc for the full callback contract, the guards, and what each column of the report is.
• ?sbc_predictive for the predictive shapes, and which one a backend should report.
• vignette("reliability-pareto-k") for the single-fit reliability band this screens against: the outer Pareto-k-hat, gamma_3, and the combined verdict.
• ?tulpa_posterior_draws for the mixture sampler the fixed-effect arms are built on, and ?tulpa_psis for the Pareto-smoothing core the k-hat uses.

References

Talts, S., Betancourt, M., Simpson, D., Vehtari, A. and Gelman, A. (2018). Validating Bayesian inference algorithms with simulation-based calibration. arXiv:1804.06788.

Sailynoja, T., Schmitt, M., Buerkner, P.-C. and Vehtari, A. (2026). Posterior SBC: simulation-based calibration checking conditional on data. Statistics and Computing 36, 78. doi:10.1007/s11222-026-10825-9.

Gneiting, T. and Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association 102, 359-378.