Martin Modrak writes:
We (Martin Modrák, Angie Moon, Shinyoung Kim, Paul Bürkner, Niko Huurre, Kateřina Faltejsková, Andrew Gelman, and Aki Vehtari) have a new preprint on SBC: Simulation-Based Calibration Checking for Bayesian Computation: The Choice of Test Quantities Shapes Sensitivity.
For those, who are unaware, simulation-based calibration checking (SBC) is a method where you use simulated datasets to verify that you implemented you model correctly and/or that your sampling algorithm work. It has been known and used for a while, but was considered to have a few shortcomings, which we try to address. See the preprint or the documentation for the SBC package for a precise description of how that works.
The gist of our preprint is that strength of SBC depends on the choice of “test quantities” for which you compute ranks. The default approach is to take all individual parameter values. This is already a very useful check, but it leaves a big class of bugs you can’t detect (e.g. when posterior is just the prior distribution). However, when you add derived test quantities, combining the parameters with the simulated data, you can (in theory) detect any problem! (Yaaaay!) But in theory you may need infinitely many quantities :-(.
In practice, it seems to be quite sufficient to add just a few additional test quantities beyond the default. In particular, our experience as well as theoretical considerations indicate that the model likelihood is very sensitive. The power of SBC is still limited by the number of simulated datasets you can reasonably fit which primarily limits how big discrepancies can go undetected.
More generally, we provide a firmly grounded theoretical analysis of SBC which will hopefully help others to build better intuition on how and why it works in practice. Notably, we make it clear that SBC does not check whether “data-averaged posterior equals the prior” as is sometimes claimed in the literature (and as I also was convinced just a couple months ago :-D )
The SBC package supports all of the ideas discussed in the preprint. I personally now use SBC for almost all my Stan work from the get go – although there is some extra work to setup SBC, you end up detecting bugs early and thus save time in the long run. If you want to see an example of simulation-driven model development workflow, check out the Small model implementation workflow.
Looking forward to your feedback on the preprint (as well as the package).
I agree that simulation-based calibration checking is an important part of statistical workflow. I often do a simpler version where I just set parameters to reasonable values, then simulate data, then fit the model to the simulated data to check to see if I can recover the parameter values to within a reasonable level of accuracy. But it’s definitely cleaner to do this by drawing from the prior.
When Cook and Rubin first suggested this idea around 2004, it was common to fit Bayesian models using hand-coded Gibbs and Metropolis algorithms, and the point of SBC was to find errors in the code that you’d write. Nowadays with Stan we’re not so worried about bugs in the sampling code, but we are concerned with possible poor mixing of the HMC and poor identification of the parameters in the model. The big idea here, though, is, as Martin writes above, to incorporate this sort of checking into workflow.
Also, for implementation you’d like to have this done in parallel, which means it would be good to have a setup where you can access 100 processors directly. I’m not sure if that’s in Martin’s workflow linked above. I guess there must be some way to link Rstudio directly with AWS so that you just put in your credit card number and it accesses multiple processors as needed?
> Nowadays with Stan we’re not so worried about bugs in the sampling code, but we are concerned with possible poor mixing of the HMC and poor identification of the parameters in the model.
We (or at least me) are also very worried that our Stan code has a bug and SBC let’s us check this.
As for computation, for most applications I work with, I get enough power from my workstation. Current implementation in the SBC package won’t play nicely with cluster for Stan models, but that’s more of a Stan problem. E.g. JAGS models are easier to transfer and thus should work withSBC and any cluster you can connect to with the future package.
Martin:
Yes, the Stan code can have bugs. That happens to me all the time. It would be a bug in the model, not in the sampling code, that’s all.
Also, currently I do most things on my laptop. But I’d like to do SBC routinely, and for that it would be convenient to be able to access a parallel server in some automatic way without requiring lots of extra programming and setup on my part.
They don’t provide you access to an HPC cluster at Columbia? I imagine that’s less expensive.
But as for AWS, you can always use AWS Sagemaker with a personal account, which has RStudio notebooks as an option. I don’t think any of their built-ins have a stan environment set up though. It’s easy enough to choose what box you want to run on, but I think most of their “scaled up” options are “many GPUs” rather than the “many CPUs” that currently work well with parallel chain MCMC.
Somebody:
I do have access to an HPC cluster at Columbia. I just don’t know how to use it. What I want is an example setup of SBC using a cluster, so I can just enter a password or credit card number or whatever and do it.
This reminds me of the famous quote in formal verification of computer programs: “Attempting to prove any nontrivial theorem about your program will expose lots of bugs: The particular choice of theorem makes little difference!”
Can someone talk me through the idea here in words? You basically run the generative model forward to sample a theta from the prior and then given that theta produce a simulated dataset, then you run the fitting procedure that samples from the posterior for theta…You do this whole procedure many times… if you do it many times and combine all the samples from the posteriors into one big mixture model, you’ll get back the prior? basically you can use your actual samples from the prior and the mixture model of all the posteriors and compare them to each other? Or am I misunderstanding?
Daniel:
Yup. The easiest place to start might be our 2006 paper: it’s not perfect and we’ve moved on, but it has a clear description of the basic idea.
I think what you describe is close in terms of overall procedure (run the generative model forward, then fit, repeat multiple times), but the actual check is different, what you described sounds like the “data-averaged posterior = prior” criterion, which seems less useful as discussed in the preprint. (also putting all the posterior samples in one bag will IMHO not do the trick, you’d need to just take a single sample from each fit).
SBC takes the rank of the prior sample within the posterior samples separately for each fit and tests that the ranks are (discrete) uniform. So in SBC you _condition_ on each specific dataset, instead of marginalizing over them (as the data-averaged posterior does). I.e. for data-averaged posterior a huge overestimate in one simulated dataset can be averaged out by underestimate in another dataset. In SBC, each will contribute an extreme low/high rank and thus not cancel out. The difference is actually quite a big deal in practice.
Does that make sense?
> when you add derived test quantities,…infinitely many quantities
This line could be alarming: having multiple test statistics or even to optimize the test statistics requires a multiple-testing framework, which is yet to be added.
Better handling of multiple testing indeed could help in a few cases, but I think that fear is overblown. The difference is that with SBC we are trying to do a severe test of the hypothesis that our model/sampling arg works. So false negatives (not discovering a bug) are much bigger issue than false positives. More test quantities can incur more false positive, but the worst case is that we end up doing some extra scrutiny of our code, not discover anything, run more simulations and see the false positive disappear.
Also in practice the test quantities are not independent, and non-uniformity in one quantity is positively correlated with non-uniformity in others, so the increase in false positives from multiple testing is quite small.
Would that convince you?