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? 
 A: 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.
A: 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.
A: 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.
