# How to use simulation to check the correctness of my Bayesian model?

To demonstrate the correctness of a frequentist estimator, it is common to simulate an experiment N times (with N being large), then show that 95% of the resulting N confidence intervals cover the true parameter values.

What's the equivalent simulation exercise for Bayesian model? The Bayesian credible interval quantifies my belief about the parameter given this one particular dataset that I got, so it doesn't make sense to simulate N experiments for N new datasets. That's where I got stuck in my thinking.

What I want to achieve specifically: I want to check whether my Stan model is implemented correctly. As recommended by the Stan manual, I generate mock data with a known DGP and fit my Stan model to it. Sometimes the resulting 95% credible interval covers the true value, sometimes not. My first reaction is to re-run this process N times and to check whether 95% of N times, my credible interval covers the true value. Is this valid?

• Unfortunately not. Consider a situation where the prior is far from the true value and the data is not strong enough to overwhelm it completely. In this case, the $\alpha$-credible interval won't have $\alpha$ coverage probability. This comes about because Bayesian statistics is about optimal updating of pre-existing information rather than about estimation without any pre-existing information being taken into account. If you make your prior match the DGP, though, you should be able to proceed as you suggest. – jbowman Jul 5 '17 at 19:42
• @jbowman when you say make the prior match the GDP, does that mean that my prior should center on the true param values of the GDP? But my prior is a distribution and thus is characterized by more than just its center. How do you mean to "match the prior with the GDP"? – Heisenberg Jul 5 '17 at 20:51
• It seems to me the two problems are the same. In the frequentist setting we have to study a range of plausible parameter values and examine the entire landscape of results, because we're not willing to stipulate a prior probability distribution. In the Bayesian setting you don't have to study a grid of parameter values, because you simply let your prior distribution choose them and average the results. To do Bayesian stats well, though, you need to evaluate the sensitivity to the prior, so maybe you should be exploring a grid of hyperparameters--and that's just like studying parameters. – whuber Jul 5 '17 at 22:17
• @whuber - I have to confess I sometimes think Bayesian statistics is what you should do when you want to update prior information and frequentist statistics what you should do when you don't want to include any prior information except for the minimum necessary to specify a likelihood function, parameter space, and perhaps a null hypothesis. I know there's a huge literature on noninformative priors, but the various methodologies generate inconsistent results and somehow don't seem as clean to me as "I'm going to do a Wilcoxon rank-sum test". Probably just lazy thinking on my part. – jbowman Jul 6 '17 at 2:20

It sounds as if you are looking for the procedure described in this paper, which is implemented in the BayesValidate R package and the pp_validate function in the rstanarm R package. Briefly, you draw repeatedly from the prior predictive distribution of the model to create simulated datasets, condition on each simulated dataset to draw from its posterior distribution, and use the quantiles of these posterior distributions to conduct a statistical test of the null hypothesis that the software is correct.

I think there's confusion, there is not really a "true value" in Bayesian inference since the parameter of inference is random. You may have to specify what form of credible interval you constructed from the posterior so we can understand what value it attempts to capture, did you use highest modal probability for the Bayesian "MLE"? did you use a symmetric interval about the posterior mean which contains 95% of the posterior? Did you choose the empirical 2.5 and 97.5 quantiles for the posterior median?

Once you specify this, you can use simulation and replicate Gibb's sampling for each simulated dataset to sample from the posterior and calculate the respective statistic represented by its credible interval.

• In this particular case, I use the empirical 2.5 and 95.5 quantile for the posterior median. However, I don't think my choice of the credible interval the issue. Are you recommending that I generate 100 datasets, fit my model 100 times and get out 100 credible intervals? Is that valid given that in Bayesian inference the data is fixed, while the parameters are random? – Heisenberg Jul 5 '17 at 19:47
• @Heisenberg no, for each simulated dataset, obtain an estimate of the posterior median. 95% of the posterior medians should be contained in the 95% credible CI you generated from the source dataset if the assumptions are met. – AdamO Jul 5 '17 at 19:50

The simulation exercise just tests whether your program has no bugs, that is, whether the algorithm you implemented correctly reflects the mathematics you want to use.

Bayesian models are not corrected: they are updated – that's an important point of the theory. There's no way to "test" the model: to do that we'd need to know what future data we'll get – which is exactly our original problem and what we're trying to guess!

There's of course the point that the models we use in our algorithms are often very simplified pictures of the complex beliefs that lie in the backs of our minds – that's why they're called "models". In this case there's an approximate update process going on in our mind which cannot be reflected in the simplified model we put on paper. In this case we start over with a new model, a new picture, which better reflects our complex beliefs. For example we use an exchangeable model for a sequence of coin tosses, and suddenly we start to see an alternating pattern in the outcomes: tail, head, tail, head, tail, head, tail, head,... The exchangeable model was too simple and will never be able to update to reflect the alternating pattern. But in our mind we're always watching for such kinds of patterns, even when we barely believe in them. So the alternating data will start increasing our belief in a non-exchangeable pattern, and we will have to cross out the one on paper, which was only a simplified picture, and replace it with a different simplified picture, for example a time-series model.

But neither the model on paper or the complex beliefs in your minds can be "tested" – they're simply updated.