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.

Installation

To install the development version of SBC, run

devtools::install_github("hyunjimoon/SBC")

Quick tour

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 generated quantities to guard against known limitations of vanilla SBC.

Paralellization

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.

library(future)
plan(multisession)

More information

The package vignettes provide additional context and examples. Notably:

  • The main vignette has more theoretical background and instructions how to integrate your own simulation code and models with SBC.
  • Small model workflow discusses how SBC integrates with model implementation workflow and how you can use SBC to safely develop complex models step-by-step.

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.

FAQ

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.

Acknowledgements

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.