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 generated 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:

- 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.

- Theoretical support
- Validating Bayesian Inference Algorithms with Simulation-Based Calibration Talts, Betancourt, Simpson, Vehtari, Gelman, 2018
- Graphical Test for Discrete Uniformity and its Applications in Goodness of Fit Evaluation and Multiple Sample Comparison Säilynoja, Bürkner, Vehtari, 2021
- Bayesian Workflow, Gelman et al., 2020
- Toward a principled Bayesian workflow in cognitive science Schad, Betancourt, Vasishth, 2021
- Bayes factor workflow Schad, Nicenboim, Bürkner, Betancourt, Vasishth, 2021

- Application support

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.