The Bayesian Workflow: Why Fitting a Posterior Is the Easy Part

In January 2017 Michael Betancourt fitted one of the smallest hierarchical models in statistics, eight SAT-coaching experiments with a population mean and a between-school standard deviation \(\tau\), and got an answer that looked fine. The 700-draw chain reported \(\hat{R} = 1.01\) for \(\tau\) and a 2.5% quantile of 0.89 (Betancourt, 2017, Diagnosing Biased Inference with Divergences). Under the same priors, direct numerical integration (done for this article) puts the 2.5% quantile at 0.12 and 17.8% of the posterior mass below 0.89. The sampler had silently cut off the region where the schools are nearly identical. Its only complaint was 12 divergent transitions, 1.7% of iterations.

That gap is why Bayesian practice now calls itself a workflow. Probabilistic programming made sampling cheap, not trustworthy.

Why this matters: A posterior is a computed object, and the computation can fail without raising an exception. The checks in this article are the difference between reporting an uncertainty estimate and reporting an artefact of the sampler. Anyone fitting models in Stan, PyMC or NumPyro needs to know what each check detects and what it misses.

TL;DR

• Bayesian workflow is iterative model building, checking, computational troubleshooting and comparison, not one fit (Gelman et al., 2020, arXiv:2011.01808).
• Independent \(N(0,1)\) priors on 15 standardised logistic-regression coefficients imply that about 45% of simulated success probabilities fall outside 0.05 to 0.95 (our calculation). Only a prior predictive simulation shows it.
• The recommended \(\hat{R}\) threshold tightened from 1.1 to 1.01 in 2021, alongside rank normalisation, folding, at least four chains and rank-normalised ESS above 400 (Vehtari et al., 2021, arXiv:1903.08008).
• Leapfrog on a Gaussian of scale \(\tau\) is stable only when the step size \(\epsilon < 2\tau\). At \(\epsilon = 0.3\) and \(\tau = 0.1\) the energy error runs 5.9, 324, 15,268 in successive steps. That is what a divergence is.
• Raising adapt_delta from 0.8 to 0.99 left the centred model's divergence count nearly constant and shrank, but never removed, its bias. The non-centred model put the 2.5% quantile at 0.11 and reached 0 divergences at 0.90.
• Simulation-based calibration makes "is my sampler right?" testable: ranks of prior draws among posterior draws are uniform on \(\{0, \dots, L\}\) exactly when inference is correct (Talts et al., 2018, arXiv:1804.06788).

At a Glance

flowchart LR
    A["Generative model"] --> B["Prior predictive check"]
    B --> C["Fit simulated data"]
    C --> D["Fit real data"]
    D --> E["Computation diagnostics"]
    E --> F["Posterior predictive check"]
    F --> G["PSIS-LOO comparison"]
    G --> H["Report or expand"]
    B -.->|"implausible data"| A
    C -.->|"SBC ranks not uniform"| A
    E -.->|"divergences, high R-hat"| A
    F -.->|"misfit"| A

    classDef blue fill:#1e40af,stroke:#3b82f6,stroke-width:1px,color:#fff
    classDef purple fill:#6d28d9,stroke:#a78bfa,stroke-width:1px,color:#fff
    classDef teal fill:#0e7490,stroke:#22d3ee,stroke-width:1px,color:#fff
    classDef amber fill:#b45309,stroke:#fbbf24,stroke-width:1px,color:#fff

    class A blue
    class B,C,E,F,G amber
    class D purple
    class H teal

Four arrows point backwards; most of the time is spent on them.

From One Chain to a Workflow

The eight-schools data come from Donald Rubin's analysis of eight randomised experiments on SAT-Verbal coaching, used to show how a hierarchical model sits between treating schools as identical and treating them as unrelated (Rubin, 1981, Journal of Educational Statistics 6(4)). The data, \(y = (28, 8, -3, 7, -1, 1, 18, 12)\) with standard errors \(\sigma = (15, 10, 16, 11, 9, 11, 10, 18)\), became the canonical stress test for hierarchical samplers.

A decade later Gelman and Rubin proposed running several chains from over-dispersed starts and comparing between-chain to within-chain variance, the potential scale reduction factor (Gelman and Rubin, 1992, Statistical Science 7(4)). The folk threshold of 1.1 stuck for almost thirty years.

timeline
    title From a single diagnostic to a Bayesian workflow
    1981 : Rubin analyses the eight schools experiments
    1992 : Gelman and Rubin propose the potential scale reduction factor
    2006-2007 : Cook, Gelman and Rubin validate software with posterior quantiles
              : Papaspiliopoulos, Roberts and Skold formalise centring versus non-centring
    2013-2014 : Betancourt and Girolami analyse HMC on hierarchical funnels
              : Hoffman and Gelman publish the No-U-Turn Sampler
    2017-2018 : PSIS-LOO makes cross-validation routine
              : Talts et al. publish simulation-based calibration
    2020-2021 : Gelman et al. write down the Bayesian workflow
              : Rank-normalised R-hat moves the threshold to 1.01
    2022-2024 : Modrak et al. add data-dependent test quantities to SBC
              : Pareto smoothed importance sampling appears in JMLR

Then the samplers changed. Hamiltonian Monte Carlo with the No-U-Turn Sampler removed most hand tuning (Hoffman and Gelman, 2014, JMLR 15) and brought a diagnostic random-walk methods never had, the divergent transition, whose natural trigger is the funnel geometry that Betancourt and Girolami analysed in hierarchical models (Betancourt and Girolami, 2013, arXiv:1312.0906).

Cook, Gelman and Rubin proposed validating software by simulating parameters and data, fitting, and asking where the truth lands in the posterior (Cook, Gelman and Rubin, 2006, JCGS 15(3)). By 2020 ten authors had assembled these pieces into a 77-page account with an uncomfortable claim: the model you report is one of many you fitted, and the discarded ones are part of the analysis.

[IMAGE: The eight schools estimates as points with plus or minus two standard error bars, school A at 28 with an interval from -2 to 58, and a horizontal line at the pooled estimate 7.7. Caption: "Every school's interval contains the pooled estimate. That is why tau is hard to pin down and the posterior is funnel-shaped."]

How the Workflow Actually Works

The model is a program that generates data

A Bayesian model is a joint distribution \(p(y, \theta) = p(y \mid \theta)\, p(\theta)\), and read generatively it is a program: draw parameters, then data. For eight schools, with the priors Betancourt's case study used:

where \(\mu\) is the average coaching effect, \(\tau\) the spread of true effects across schools, \(\theta_j\) school \(j\)'s true effect and \(\sigma_j\) its known standard error.

Prior predictive checks: priors act on data

A prior predictive check runs the program forward with no observed data and asks whether the simulated datasets look possible. It exists because priors combine invisibly. Gelman et al. show independent \(N(0,1)\) logistic-regression coefficients pushing the simulated mean outcome toward 0 or 1 as the number of binary covariates grows from 2 to 4 to 15 (Gelman et al., 2020, Section 2.4).

The mechanism is one line (our illustration, with standardised rather than binary covariates). With \(p\) covariates and \(\beta_k \sim N(0,1)\), the linear predictor \(\eta = \sum_k \beta_k x_k\) has prior variance

A probability leaves \([0.05, 0.95]\) when \(|\eta| > \operatorname{logit}(0.95) = 2.944\). With \(p = 2\), \(\text{sd}(\eta) = 1.41\) and \(P(|Z| > 2.08) \approx 0.04\). With \(p = 15\), \(\text{sd}(\eta) = 3.87\) and \(P(|Z| > 0.76) \approx 0.45\). The same "weak" prior now asserts that nearly half of all observations are near-certain. Fitting would never reveal it.

What a divergence is, derived

HMC treats \(U(\theta) = -\log p(\theta \mid y)\) as potential energy, adds momentum \(r\), and simulates \(H = U + \tfrac{1}{2}r^\top r\) with the leapfrog integrator. Exact flow conserves \(H\); leapfrog conserves it approximately, and the approximation fails in an analysable way.

Take one school effect in the funnel with \(\mu\) fixed: \(U(\theta) = \theta^2/(2\tau^2)\), an oscillator with frequency \(\omega = 1/\tau\). One leapfrog step of size \(\epsilon\) is the linear map

with determinant 1 (leapfrog preserves volume) and trace \(2 - \epsilon^2\omega^2\). The eigenvalues solve \(\lambda^2 - (2 - \epsilon^2\omega^2)\lambda + 1 = 0\). While the trace lies in \((-2, 2)\) they are complex with modulus 1 and energy error stays bounded. Once the trace drops below \(-2\), one eigenvalue exceeds 1 in modulus and the error grows geometrically. Stability requires

Stan adapts one step size during warmup to a target acceptance rate (adapt_delta, default 0.8), tuned to the bulk where \(\tau\) is several units. In the funnel's neck \(\epsilon\omega\) crosses 2 and the simulated energy explodes. Stan flags a divergence when \(H\) departs from its starting value by more than \(10^3\) (Stan Reference Manual, MCMC Sampling), stops that trajectory and carries on. Nothing crashes. The chain simply never samples the region it could not integrate, so divergences imply bias, not just inefficiency.

Funnel geometry and reparameterisation

Stan samples \(\nu = \log\tau\), adding a Jacobian term \(\nu\). The centred model's log density in \((\theta, \nu)\) is

The curvature along each \(\theta_j\) is \(1/\tau^2\), which across the exact 95% interval for \(\tau\) (0.12 to 11.93) varies by a factor of nearly \(10^4\). And with every \(\theta_j\) near \(\mu\), the term \(-(J-1)\nu\) grows without bound as \(\nu \to -\infty\): for \(J = 8\) density rises like \(e^{7|\nu|}\) into vanishing volume. That neck holds real mass (10.1% of the posterior lies below \(\tau = 0.5\)) and is exactly where stability fails.

The non-centred parameterisation samples \(\tilde\theta_j \sim N(0,1)\) and sets \(\theta_j = \mu + \tau\tilde\theta_j\). Prior curvature along every \(\tilde\theta_j\) is now 1 whatever \(\tau\) is; the likelihood couples them with curvature \(\tau^2/\sigma_j^2\), small exactly when data are weak. Same model, different coordinates.

Papaspiliopoulos, Roberts and Sköld showed centred parameterisations mix well when each group's data are informative relative to the population spread, and non-centred ones when data are weak (Papaspiliopoulos, Roberts and Sköld, 2007, Statistical Science 22(1)). Eight schools, with standard errors of 9 to 18 against a \(\tau\) of a few points, is firmly non-centred territory. Thousands of students per school would reverse the advice.

[IMAGE: Two scatter plots of posterior draws. Left, centred: theta_1 against log tau forming a funnel, green divergent draws clustered at the neck, empty region below log tau of minus 1. Right, non-centred: tilde theta_1 against log tau, a rectangular cloud reaching log tau of minus 5. Caption: "The funnel belongs to the parameterisation, not the model."]

R-hat, rebuilt

Split-\(\hat{R}\) halves each chain into \(M\) sequences of \(N\) draws. With within-sequence variance \(W\) and \(B = \frac{N}{M-1}\sum_m (\bar\theta_m - \bar\theta)^2\),

Vehtari et al. found it misses chains with equal means but different scales and misbehaves for heavy tails; in their simulations it failed to flag chain sets built to be far from mixed, such as one chain with a third of the others' variance (Vehtari et al., 2021, Bayesian Analysis 16(2)).

The fixes transform draws before the same formula. Rank normalisation replaces each of the \(S\) pooled draws by its rank \(r\) and maps it to \(z = \Phi^{-1}\big((r - 3/8)/(S + 1/4)\big)\), defined for any distribution. Folding applies the statistic to \(|\theta - \operatorname{median}(\theta)|\), turning a scale difference into a location difference. Report the maximum of the two, use the draws only if \(\hat{R} < 1.01\), run at least four chains, and require rank-normalised ESS above 400.

Effective sample size is \(\text{ESS} = S / (1 + 2\sum_{t \ge 1}\rho_t)\) for lag-\(t\) autocorrelations \(\rho_t\), and the Monte Carlo standard error of a mean is \(\text{sd}/\sqrt{\text{ESS}}\). Bulk ESS governs means; tail ESS governs the 5% and 95% quantiles, where funnel bias lives. For the centred eight-schools model under conservative settings, Vehtari et al. report bulk and tail ESS for \(\tau\) of 67 and 82.

Simulation-based calibration: testing the computation

Bayesian inference is calibrated only on average over the prior, so SBC makes that average the test. Draw \(\tilde\theta \sim p(\theta)\), simulate \(\tilde y\), take \(L\) posterior draws, and compute

for a scalar test quantity \(f\). With exact inference, \((\tilde\theta, \theta_1, \dots, \theta_L)\) are exchangeable, so the rank is uniform on \(\{0, \dots, L\}\) (Talts et al., 2018). Repeat \(N\) times and compare each histogram bin with the 0.5% and 99.5% quantiles of \(\text{Binomial}(N, 1/(L+1))\).

The shape names the failure. \(\cup\)-shaped: posteriors too narrow. \(\cap\)-shaped: too wide. Tilted: biased. Spikes at both ends: autocorrelated draws needing thinning. Talts et al. ran it with \(L = 100\) on these very models: the centred eight-schools fit gave a clearly non-uniform histogram for \(\tau\), the thinned non-centred fit a uniform one, and Stan's ADVI a strongly biased one for the slope of a simple linear regression.

Modrák et al. showed the choice of \(f\) decides what SBC sees. Parameter-only ranks miss a posterior that collapses to the prior, which is still calibrated on average; adding the joint log likelihood as a test quantity catches it (Modrák et al., 2025, arXiv:2211.02383).

Posterior predictive checks and PSIS-LOO

Once computation is trusted, the question is adequacy. A posterior predictive check draws \(y^{\text{rep}} \sim p(y^{\text{rep}} \mid y)\) and compares \(T(y^{\text{rep}})\) with \(T(y)\). Model comparison uses leave-one-out predictive density,

where raw ratios \(r_i^s = 1/p(y_i \mid \theta^s)\) reweight full-data draws. Those ratios are heavy-tailed whenever observation \(i\) is influential, so Pareto smoothed importance sampling fits a generalised Pareto distribution to the largest and replaces them with smoothed values (Vehtari, Gelman and Gabry, 2017, Statistics and Computing 27). The fitted shape \(\hat{k}\) is itself a diagnostic, since a shape-\(k\) tail has finite moments only below order \(1/k\). Estimates are reliable when \(\hat{k} < \min(1 - 1/\log_{10} S,\ 0.7)\), and \(\hat{k} \ge 0.7\) means large bias (loo documentation; Vehtari et al., 2024, JMLR 25).

Seeing It in Motion

One transition that diverges

sequenceDiagram
    participant S as NUTS sampler
    participant L as Leapfrog integrator
    participant G as Gradient
    participant D as Divergence monitor
    S->>L: Start near log tau of -2.3 with fresh momentum
    L->>G: Request gradient
    G-->>L: Curvature about 100
    L->>D: Energy error 5.9 after step 2
    L->>D: Energy error 324 after step 3
    L->>D: Energy error 15268 after step 4
    D-->>S: Above 1000, flag divergent
    Note over S: Tree building stops, draw taken from partial trajectory
    S->>S: Record divergence and continue

The numbers come from the worked example; the control flow is Stan's.

Triage when the sampler complains

stateDiagram-v2
    [*] --> Fitted
    Fitted --> CheckDivergences
    CheckDivergences --> Clustered: concentrated in one region
    CheckDivergences --> Scattered: rare and spread out
    CheckDivergences --> CheckRhat: none
    Scattered --> RaiseAdaptDelta
    RaiseAdaptDelta --> CheckRhat: divergences vanish
    RaiseAdaptDelta --> Clustered: divergences persist
    Clustered --> Reparameterise
    Reparameterise --> Fitted
    CheckRhat --> ReviseModel: R-hat above 1.01 or ESS below 400
    CheckRhat --> PredictiveChecks: all pass
    ReviseModel --> Fitted
    PredictiveChecks --> [*]

    classDef amber fill:#b45309,stroke:#fbbf24,stroke-width:1px,color:#fff
    classDef teal fill:#0e7490,stroke:#22d3ee,stroke-width:1px,color:#fff
    class Clustered amber
    class Reparameterise amber
    class ReviseModel amber
    class PredictiveChecks teal

A smaller step fixes only scattered divergences, often false positives; clustered ones mean geometry. The Stan team notes that raising adapt_delta has become the reflexive first move and that going beyond 0.99 is seldom useful (Stan Development Team, Runtime warnings and convergence problems).

[IMAGE: Four SBC rank histograms with grey 99% binomial bands: uniform (correct), U-shaped (posterior too narrow), inverted U (too wide), ranks piled high (draws biased low). Caption: "The shape of the deviation names the failure."]

By the Numbers

Sources: Betancourt (2017) (single chains, priors as above); Vehtari et al. (2021); Stan guide; Stan Reference Manual; Talts et al. (2018); loo documentation. "Exact" rows are this article's grid over \(\tau\) with \(\mu\) integrated analytically, close to the long non-centred run. Single chains illustrate the failure; they are not benchmarks. [IMAGE: Running MCMC mean of log tau over 10,000 iterations. The centred chain oscillates, sticks, then settles visibly above a dashed reference line from a long non-centred run; the non-centred chain settles onto it quickly. Caption: "The centred chain converges, confidently, to the wrong number."]

A Concrete Example

Step 1: No pooling and complete pooling. Alone, school A estimates \(28 \pm 15\). Pooled, the precision weights \(1/\sigma_j^2\) are

with \(\sum_j y_j/\sigma_j^2 = 0.4635\), so the pooled mean is \(0.4635/0.06031 = 7.69\) with standard error \(1/\sqrt{0.06031} = 4.07\).

Step 2: Heterogeneity. \(Q = \sum_j (y_j - 7.69)^2/\sigma_j^2 = 4.71\) on 7 degrees of freedom, below its null expectation of 7 (p = 0.70). The moment estimate of \(\tau^2\) is negative and truncates to zero, yet the data are also compatible with \(\tau = 10\). That ambiguity is the funnel.

Step 3: Partial pooling given \(\tau\). School \(j\) shrinks toward \(\mu\) by \(B_j = \sigma_j^2/(\sigma_j^2 + \tau^2)\):

For school A (\(\sigma_A^2 = 225\)), with flat priors on \(\mu\) and \(\tau\) for this illustration:

Averaged over the posterior of \(\tau\) under the same flat priors (grid calculation), school A's posterior mean is about 11.4, down from 28; school C rises from \(-3\) to about 6.1. Every school's answer depends on the mass at small \(\tau\), the region the centred sampler loses.

Step 4: Why it loses it. Take one centred coordinate at \(\tau = 0.1\) (\(\omega = 10\)) with \(\epsilon = 0.3\), a step comfortable at \(\tau = 5\). Start at \(\theta = 0.1\), \(r = 1\), so \(H_0 = 0.5 + 0.5 = 1\). Now \(\epsilon\omega = 3 > 2\); the eigenvalues solve \(\lambda^2 + 7\lambda + 1 = 0\), giving \(-6.85\) and \(-0.146\), so error multiplies by about 6.85 per step with alternating sign:

Step 4 crosses 1,000 and the transition is flagged. The same recursion at \(\tau = 5\) keeps the error below 0.001 for ten steps, and at \(\tau = 1\) below 0.02.

Step 5: The cost. Under the case-study priors, 17.8% of the exact posterior lies below \(\tau = 0.89\), yet the centred chain reported 0.89 as its 2.5% quantile. About 15 percentage points of "nearly identical schools" went missing, understating shrinkage and overstating heterogeneity. The non-centred chain reported 0.11, and \(\hat{R}\) could not tell the two runs apart.

[IMAGE: Leapfrog trajectories in the theta-momentum plane at step size 0.3 for tau of 5 (smooth ellipse), 1 (slightly distorted orbit) and 0.1 (points spiralling outward with alternating sign). Caption: "A divergence is an integrator crossing epsilon times omega equals 2."]

Where It Breaks

Diagnostics detect non-convergence; none certify convergence

\(\hat{R}\), ESS and divergences flag only problems the chains encountered. Chains all started in one mode of a multimodal posterior pass every check. In the workflow paper's planetary-motion example chains stuck in different modes, and the fix was understanding the likelihood surface, not running longer (Gelman et al., 2020, Section 11).

The threshold debate is not settled

Vats and Knudson argue \(\hat{R}\) is a one-to-one function of effective sample size once variance is estimated properly, so a threshold should follow from the precision you need rather than be fixed (Vats and Knudson, 2021, Statistical Science 36(4)). Both camps agree \(\hat{R}\) alone is insufficient; they disagree about whether a universal cutoff means anything.

Raising adapt_delta hides bias

A smaller step lets trajectories go deeper into the neck. In the case study the divergence count stayed nearly constant from 0.80 to 0.99 while \(\tau\)'s distribution moved toward the truth without reaching it. Where the count does fall, a run down to a handful of divergences looks fixed and can still be biased.

Non-centring can create the problem it solves

With data-rich groups the likelihood forces \(\mu + \tau\tilde\theta_j \approx y_j\), a narrow ridge in non-centred coordinates, so mixed datasets need some groups centred and others not.

SBC is expensive and prior-sensitive

SBC needs \(N\) full fits. Weakly informative priors also generate extreme datasets that break samplers which handle realistic data; the workflow paper concedes SBC "clashes somewhat" with how practitioners set priors.

PSIS-LOO answers a predictive question

Gronau and Wagenmakers built cases where the simpler model is exactly true and showed LOO's support for it stays bounded even with infinite data (Gronau and Wagenmakers, 2019, Computational Brain & Behavior 2). Vehtari, Simpson, Yao and Gelman replied that bounded support is expected of a predictive criterion, that true-model-in-the-set cases are rare, and that LOO is not a hypothesis test (Vehtari et al., 2018, arXiv:1810.05374). Mechanically, with one parameter per observation, removing observation \(i\) reshapes the posterior, pushes \(\hat{k}\) above 0.7, and forces refitting or K-fold cross-validation.

The folk theorem cuts both ways

"When you have computational problems, often there's a problem with your model" (Gelman et al., 2020, Section 5.1) is the right first hypothesis, not a guarantee. Bending a justified model to please a sampler trades computational bias for modelling bias.

Alternative Designs

Neither of the first two rows is a default; the ratio of per-group data precision to population spread decides. Marginalisation is underused: for Gaussian eight schools it gives this article's exact reference.

How It Is Used in Practice

Stan, PyMC and NumPyro run NUTS-descended HMC; posterior, loo and ArviZ report rank-normalised \(\hat{R}\), bulk and tail ESS, divergences and Pareto \(\hat{k}\).

Before real data. Plot prior predictive summaries on the outcome scale and confirm recovery on simulated datasets. For a model refitted weekly, run SBC once with the log likelihood among the test quantities; the cost amortises over every future fit.

On every fit. Four or more chains from dispersed starts; an explanation for every post-warmup divergence; \(\hat{R} < 1.01\) and bulk and tail ESS above 400 for every reported quantity, derived ones included. Divergences clustered at small hierarchical scales point at a centred parameterisation.

After computation is trusted. Choose predictive statistics for the decision, such as tail frequencies for a risk model. Report \(\text{elpd}_{\text{loo}}\) differences with standard errors, refit observations above the \(\hat{k}\) bound, and keep discarded models as evidence. SBC is the expensive step teams skip: defensible for a familiar model class, indefensible for a new one.

[IMAGE: Checklist-style flow for one production refit: prior predictive plot, four chains, divergence scatter against log tau, R-hat and ESS table, posterior predictive tail check, LOO comparison with k-hat bars. Caption: "The routine costs minutes; the SBC run behind it is done once per model."]

Insights Worth Remembering

1. A clean \(\hat{R}\) is compatible with a badly biased posterior. The centred chain read 1 over 10,000 draws while losing about 15 percentage points of mass. Divergences and a collapsing ESS per iteration were the signals.

2. Divergences are a stability limit, not a bug. Leapfrog in curvature \(1/\tau^2\) is stable only for \(\epsilon < 2\tau\). With that inequality, the funnel, step-size adaptation and non-centring become one idea.

3. Weak priors on parameters can be strong priors on data. Fifteen \(N(0,1)\) coefficients put about 45% of simulated probabilities outside 0.05 to 0.95.

4. Parameterisation is a large computational lever that leaves the model unchanged. Over 10,000 draws, centred eight schools gave ESS 376 for \(\tau\) and a biased tail; non-centred gave 6,083 and the right quantiles.

5. More iterations do not fix geometry. From 700 to 10,000 centred draws the 2.5% quantile moved from 0.89 to 0.67, against 0.12.

6. SBC tests the algorithm and predictive checks test the model. A calibrated sampler for a misspecified model passes SBC; a biased sampler can pass a predictive check.

7. Thresholds are conventions, not theorems. 1.01, 400 and 0.7 came from simulation and experience. The real criterion is a Monte Carlo standard error small enough for the decision.

Open Questions

What replaces \(\hat{R}\) for many short chains? GPU samplers run hundreds of short chains, where classic \(\hat{R}\) needs long ones. Nested \(\hat{R}\) groups chains into superchains and has theoretical backing (Margossian et al., Bayesian Analysis, 2024, arXiv:2110.13017); its behaviour on large production models is still being measured, and no consensus threshold exists.

Can reparameterisation be chosen automatically? Partially non-centred forms are robust to data informativeness (Papaspiliopoulos et al., 2007); whether languages can pick them by default without new failure modes is open.

How few SBC replications are enough? The workflow paper calls this open and speculates that designed prior draws would beat random ones. Data-dependent test quantities demonstrably raise sensitivity; no principled rule for \(N\) exists.

How should workflow be recorded? Gelman et al. argue the sequence of fitted models is evidence, yet no reporting standard for discarded models and failed checks is widely adopted; whether one emerges is speculation.

Sources and Further Reading

1. Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., et al. (2020). "Bayesian Workflow." arXiv:2011.01808

2. Rubin, D. B. (1981). "Estimation in Parallel Randomized Experiments." Journal of Educational Statistics, 6(4), 377-401. doi:10.3102/10769986006004377

3. Gelman, A., & Rubin, D. B. (1992). "Inference from Iterative Simulation Using Multiple Sequences." Statistical Science, 7(4). doi:10.1214/ss/1177011136

4. Cook, S. R., Gelman, A., & Rubin, D. B. (2006). "Validation of Software for Bayesian Models Using Posterior Quantiles." Journal of Computational and Graphical Statistics, 15(3), 675-692. doi:10.1198/106186006X136976

5. Papaspiliopoulos, O., Roberts, G. O., & Sköld, M. (2007). "A General Framework for the Parametrization of Hierarchical Models." Statistical Science, 22(1), 59-73. doi:10.1214/088342307000000014

6. Betancourt, M., & Girolami, M. (2013). "Hamiltonian Monte Carlo for Hierarchical Models." arXiv:1312.0906

7. Hoffman, M. D., & Gelman, A. (2014). "The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo." JMLR, 15(47), 1593-1623. JMLR

8. Betancourt, M. (2017). "Diagnosing Biased Inference with Divergences." Stan case study. mc-stan.org

9. Vehtari, A., Gelman, A., & Gabry, J. (2017). "Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC." Statistics and Computing, 27(5), 1413-1432. arXiv:1507.04544

10. Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A. (2018). "Validating Bayesian Inference Algorithms with Simulation-Based Calibration." arXiv:1804.06788

11. Gronau, Q. F., & Wagenmakers, E.-J. (2019). "Limitations of Bayesian Leave-One-Out Cross-Validation for Model Selection." Computational Brain & Behavior, 2(1), 1-11. doi:10.1007/s42113-018-0011-7. Reply: Vehtari, A., Simpson, D. P., Yao, Y., & Gelman, A. arXiv:1810.05374

12. Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). "Rank-Normalization, Folding, and Localization: An Improved R-hat for Assessing Convergence of MCMC." Bayesian Analysis, 16(2). doi:10.1214/20-BA1221; arXiv:1903.08008

13. Vats, D., & Knudson, C. (2021). "Revisiting the Gelman-Rubin Diagnostic." Statistical Science, 36(4). doi:10.1214/20-STS812

14. Vehtari, A., Simpson, D., Gelman, A., Yao, Y., & Gabry, J. (2024). "Pareto Smoothed Importance Sampling." JMLR, 25(72), 1-58. arXiv:1507.02646

15. Modrák, M., Moon, A. H., Kim, S., Bürkner, P., et al. (2025). "Simulation-Based Calibration Checking for Bayesian Computation: The Choice of Test Quantities Shapes Sensitivity." Bayesian Analysis, 20(2). doi:10.1214/23-BA1404; arXiv:2211.02383

16. Stan Development Team. "Runtime warnings and convergence problems." mc-stan.org

17. Stan Development Team. Stan Reference Manual, "MCMC Sampling." mc-stan.org

18. Stan Development Team. "Diagnostics for Pareto smoothed importance sampling." loo package documentation. mc-stan.org/loo

Free to read, no ads, no sign-up. If it was useful you can buy me a coffee.