SBC provides tools to validate your Bayesian model and/or a sampling algorithm via the self-recovering property of Bayesian models. This package lets you run SBC easily and perform postprocessing and visualisations of the results to assess computational faithfulness.
To install the development version of SBC, run
devtools::install_github("hyunjimoon/SBC")To use SBC, you need a piece of code that generates simulated data that should match your model (a generator) and a statistical model + algorithm + algorithm parameters that can fit the model to data (a backend). SBC then lets you discover when the backend and generator don’t encode the same data generating process (up to certain limitations).
For a quick example, we’ll use a simple generator producing normally-distributed data (basically y <- rnorm(N, mu, sigma)) with a backend in Stan that mismatches the generator by wrongly assuming Stan parametrizes the normal distribution via precision (i.e. it has y ~ normal(mu, 1 / sigma ^ 2)).
library(SBC)
gen <- SBC_example_generator("normal")
# interface = "cmdstanr" or "rjags" is also supported
backend_bad <- SBC_example_backend("normal_bad", interface = "rstan")Note: Using the cmdstanr interface, a small number of rejected steps will be reported. Those are false positives and do not threaten validity (they happen during warmup). This is a result of difficulties in parsing the output of cmdstanr. We are working on a resolution.
You can use SBC_print_example_model("normal_bad") to inspect the model used.
We generate 50 simulated datasets and perform SBC:
ds <- generate_datasets(gen, n_sims = 50)
results_bad <- compute_SBC(ds, backend_bad)The results then give us diagnostic plots that immediately show a problem: the distribution of SBC ranks is not uniform as witnessed by both the rank histogram and the difference between sample ECDF and the expected deviations from theoretical CDF.
plot_rank_hist(results_bad)
plot_ecdf_diff(results_bad)We can then run SBC with a backend that uses the correct parametrization (i.e. with y ~ normal(mu, sigma)):
backend_sd <- SBC_example_backend("normal_sd", interface = "rstan")
results_sd <- compute_SBC(ds, backend_sd)
plot_rank_hist(results_sd)
plot_ecdf_diff(results_sd)The diagnostic plots show no problems in this case. As with any other software test, we can observe clear failures, but absence of failures does not imply correctness. We can however make the SBC check more thorough by using a lot of simulations and including suitable derived quantities to guard against known limitations of vanilla SBC.
The examples above are very fast to compute, but in real use cases, you almost certainly want to let the computation run in parallel via the future package.
The package vignettes provide additional context and examples. Notably:
Currently SBC supports cmdstanr, rstan, and brms models out of the box. With a little additional work, you can integrate SBC with any exact or approximate fitting method as shown in the Implementing backends vignette.
Developing and maintaining open source software is an important yet often underappreciated contribution to scientific progress. Thus, whenever you are using open source software, please make sure to cite it appropriately so that developers get credit for their work.
When using SBC, please cite the following publication:
Modrák, M., Moon, A. H., Kim, S., Bürkner, P., Huurre, N., Faltejsková, K., … & Vehtari, A. (2023). Simulation-based calibration checking for Bayesian computation: The choice of test quantities shapes sensitivity. Bayesian Analysis, advance publication, DOI: 10.1214/23-BA1404.
Further, SBC relies on several other R packages and, of course, on R itself. To find out how to cite R and its packages, use citation function. SBC specifically rely on posterior for manipulating posterior draws and future for parallel processing. Also, do not forget to cite the probabilistic inference tool you use as backend (e.g. Stan, JAGS, brms, …)
How does calibration relate to prediction accuracy?
Comparing the ground truth and the simulated result is a backbone of calibration and comparison target greatly affects the calibrated (i.e. trained) result, similar to reward in reinforcement learning. In this sense, if the U(a(y), theta) term is designed for prediction, the model will be calibrated to have best predictive result as possible.
Development of this package was supported by ELIXIR CZ research infrastructure project (Ministry of Youth, Education and Sports of the Czech Republic, Grant No: LM2018131) including access to computing and storage facilities.