Shinyoung Kim, Angie Hyunji Moon, Martin Modrák, and Teemu Säilynoja prepared this R package for simulation-based calibration (SBC)! They write:

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 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 variance (i.e. it has y ~ normal(mu, sigma ^ 2)).

I decided to try it out. The first step was to open R and install the devtools package. From there I copied the code from the above-linked webpage, just using cmdstanr instead of rstan:

devtools::install_github("hyunjimoon/SBC")
library(SBC)
gen <- SBC_example_generator("normal")
backend_var <- SBC_example_backend("normal_var", interface = "cmdstanr")
SBC_print_example_model("normal_var")
ds <- generate_datasets(gen, n_sims = 50)
results_var <- compute_SBC(ds, backend_var)

It gives some output:

I then continue with the code from the instructions page:

plot_rank_hist(results_var)
plot_ecdf_diff(results_var)

Which gives these two pairs of plots:

I then continued by running SBC with their backend that uses the correct parametrization (i.e. with y ~ normal(mu, sigma)):

backend_sd <- SBC_example_backend("normal_sd", interface = "cmdstanr")
results_sd <- compute_SBC(ds, backend_sd)

And this is what comes out. The result is not completely clear:

So let's do the graphs:

plot_rank_hist(results_sd)
plot_ecdf_diff(results_sd)

Kim et al. write:

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.

More information is in the vignettes for the package. Martin also posted this description in the Stan forums. He describes it as "the closest we have to unit tests for statistical models."

And, as usual, everything is free and transparent so this code should be easily adaptable to other probabilistic programming languages.

This is great, and we've really gone far since Samantha Cook first had this idea as part of her Ph.D. thesis, which back in 2006 became this article and this R package.

I'm just kinda bothered that it spits out a warning even for the model that works fine. I wonder if they could fix that?