Fitting and Diagnosing Stan Models

Follows Bayesian Workflow Part 3 (ch. 11–13). Two rules organize everything here:

• Fit fast, fail fast. Run short chains on a new model, look at the diagnostics, and fix the model before spending compute on precise inference for something you may discard.
• The folk theorem of statistical computing. When you have computational problems, there is usually a problem with your model. The first instinct should be to look at the posterior and the model, not to raise adapt_delta, max_treedepth, or the iteration count.

For writing the program see stan-modeling; for interface details in R see r-stan-interfaces.

Default run settings

fit = model.sample(
    data=data, chains=4, parallel_chains=4,
    iter_warmup=1000, iter_sampling=1000,
    seed=20260914, inits=0.1, adapt_delta=0.8, show_progress=False,
)
fit.save_csvfiles("results/model_v1")   # save immediately

fit <- mod$sample(data = data, chains = 4, parallel_chains = 4,
                  iter_warmup = 1000, iter_sampling = 1000,
                  seed = 20260914, init = 0.1, refresh = 500)
fit$save_object("results/model_v1.rds")

Stan's default init range, uniform(−2, 2) on the unconstrained scale, assumes unit-scale parameters. The book now recommends init = 0.1; with predictors far from unit scale even that can overflow (exp(alpha + beta * x) with x ~ 2000). The durable fix is to put data and parameters on unit scale.

Reading the diagnostics

Check in this order after every fit. Both interfaces print warnings; also call fit.diagnose() (CmdStanPy) or fit$diagnostic_summary() / fit$cmdstan_diagnose() (cmdstanr).

print(fit.diagnose())
summ = fit.summary()                     # Mean, MCSE, StdDev, 5%, 50%, 95%, ESS_bulk, ESS_tail, R_hat
print(summ.loc[summ["R_hat"] > 1.01])
print(fit.divergences, fit.max_treedepths)

fit$diagnostic_summary()                 # divergences, treedepth, E-BFMI per chain
fit$summary(NULL, posterior::default_convergence_measures())

Rank plots (az.plot_rank, mcmc_rank_overlay) are more sensitive than trace plots. Trace plots are for understanding a failure, not for certifying success.

How many iterations and digits

Decide the digits first (BW §11.4–11.6). Two significant digits is usually enough; one is fine while exploring.

1. Run with defaults; check R-hat and ESS for all quantities.
2. Confirm ESS is high enough for the diagnostics themselves to be reliable (≥ 100 per chain).
3. Look at the posterior sd of each quantity of interest and pick the digits worth reporting.
4. Check mcse_mean, mcse_q5, mcse_q95 (posterior R package: summarise_draws(draws, default_mcse_measures()); ArviZ: az.summary(..., kind="diagnostics")). Twice the MCSE is the likely variation across seeds.
5. If the last reported digit moves, report fewer digits or run more. Halving MCSE needs 4× the draws.

Rules of thumb: 100 independent draws give one significant digit for a mean; 2000 give two. Stan's NUTS often yields ESS > S/2, so 4 × 1000 draws is close to two digits for means and one to two for 5%/95% quantiles. Indicators (Pr(beta > 0)) need many more draws for the same digits. For quantities with infinite variance use median and MAD; quantile MCSE is always defined.

Failure modes and the step forward

Full worked examples in references/failure_modes.md.

Divergence rule (BW §12.3): if more than about 1% of transitions diverge, raising adapt_delta will not help. Look at pairs plots colored by divergence (az.plot_pair(dt, visuals={"divergence": True}), mcmc_pairs(np = nuts_params(fit)), mcmc_scatter(..., np = ...)) and find where they cluster.

Reparameterization toolbox

• Non-centered hierarchical: mu = mu0 + sigma0 * z; z ~ std_normal(); when groups have little data. Centered when groups have lots of data. Mixed cases have no automatic fix; parameterize groups differently if needed.
• <offset=, multiplier=> on declarations lets Stan sample on the scaled space while you write the centered model: vector<offset=mu0, multiplier=sigma0>[K] mu;
• Center and scale predictors in transformed data; back-transform coefficients in generated quantities.
• QR decomposition for correlated predictors: Q = qr_thin_Q(X) * sqrt(N-1), R = qr_thin_R(X) / sqrt(N-1), model on theta, recover beta = R \ theta.
• Well-identified quantities as the basis: start time + log duration; total variance + simplex split across variance components; for a sigmoid, the value at a reference x and the inflection point instead of raw (a, b).
• Marginalize the local parameters when possible (normal-normal, GP with normal likelihood, beta-binomial instead of logit-normal binomial); sample the hyperparameters, then draw locals in generated quantities.
• Zero-avoiding priors (lognormal, inverse-gamma) on group scales remove the funnel neck; legitimate when there is such information, but say so and run a sensitivity check if used only for speed.

See references/reparameterization.md for Stan code.

Debugging strategy

• Meet in the middle: simplify the failing model until it works; build up the working simple model until it breaks. The bug lives between.
• Simulate from the model and fit the simulated data; misspecified models are slow and simulated data separate computation from modeling problems.
• Modulate the prior: add normal(0, 100) to every parameter as a diagnostic (only matters if the posterior is improper); when freeing a fixed parameter, go normal(4, 0.1) then normal(4, 1), not straight to flat.
• Fit a subset, 200 iterations, fewer chains, while diagnosing.
• Optimize first: model.optimize() / mod$optimize() smoke-tests the program in seconds.
• Print intermediates for numerical problems: evaluate log_prob at a problematic draw (model.log_prob(params, data), fit$log_prob()), look for Inf, 0, NaN.
• Slow model checklist (BW §12.4): simulate, build up one component at a time, 200 iterations, informative priors on coefficients and group scales, question the model (interactions?), subset the data.

Approximate algorithms

Approximate algorithm = exact algorithm for an approximate model. Validate any approximation the same way as a model: simulate data with known parameters, run the pipeline, compare.

Parallelism

• Between chains: parallel_chains=4.
• Within chain: reduce_sum over the likelihood with cpp_options={"STAN_THREADS": True} and threads_per_chain=. Worthwhile when the likelihood dominates cost (large N, ODEs). See references/performance.md.
• GPU/OpenCL for large matrix and GLM ops via STAN_OPENCL.

Replicability

Bit-level reproducibility across machines is not achievable; aim for agreement within the reported digits, justified by MCSE. Fix seeds while iterating and in demos; for final results either report the seed or show that the digits are stable across seeds.

References

• references/failure_modes.md: each failure mode with the Stan code, the warning text, what the diagnostics show, and the fix.
• references/reparameterization.md: centered vs non-centered, offset/multiplier, QR, marginalization, mixture identifiability.
• references/approximate_inference.md: optimize, Laplace, Pathfinder, ADVI, divide-and-conquer, PSIS updating, model simplifications.
• references/performance.md: vectorization, GLM functions, reduce_sum, profiling, compile flags.