The core of Modern Statistical Workflow is to always simulate fake data from your generative model before even touching real data. Whenever a friend asks for help with fitting a model in Stan, my first question is: “have you estimated the model on fake data simulated from the generative process you are proposing?” More often than not, the answer is “no”, and in the process of simulating the fake data they catch their own errors.

When simulating this fake data, our parameters should be drawn from the prior. Doing so tells us when our priors are stupid. Whenever you’re having to try to defend yourself against someone who’s telling you to use some crazily diffuse prior like \(\beta \sim \mbox{Normal}(0, 1000)\) tell them to go through this exercise (we do this below). There’s also sense in drawing your parameters from the prior and then estimating your model on the resulting fake data: when your prior places little weight on the true value of a parameter, your inference (and sampling) will be compromised.

So how can we make this process easier?

The problem with this approach is that typically we have to take care of two bits of code: our simulation model and the estimation model. If we make a change in one, that change needs to be reflected in the other. That might not be an issue when your model is simple. But if you’re in the world of fitting models in Stan, then you’re trying to fit complex models. And making changes to two sets of code is a royla PITA.

Making use of Stan’s generated quantities block to simulate fake data

The really simple fix to this is to handle all modeling, both simulation and estimation, in the one piece of code. We do this by only evaluating the likelihood if we tell it to. If a likelihood is not evaluated, then our posterior is equal to the prior. In this case, our posterior predictive draws are actually prior predictive draws—precisely the fake data we want.

As a very simple example, let’s put together a tiny classification model. Say that our outcome \(y\) is distributed according to a categorical distribution with probabilities vector \(\mbox{softmax}(\theta)\)

\[ y_{i} \sim \mbox{Categorical}(\mbox{softmax}(\theta)) \]

To identify the model, we’ll set the first element of \(\theta\) to 0. Over the rest of the parameters we’ll give a very diffuse prior, say, \(\theta_{-1} \sim \mbox{normal}(0, 100)\).

We code this up in Stan as:

// saved as categorical_model.stan
data {
  int N;
  int P; // number of categories to be estimated
  int y[N]; // outcomes
  int<lower = 0, upper = 1> run_estimation; // a switch to evaluate the likelihood
  real<lower = 0> prior_sd; // standard deviation of the prior on theta
parameters {
  vector[P-1] theta_raw;
transformed parameters {
  vector[P] theta;
  theta[1] = 0.0;
  theta[2:P] = theta_raw;
model {
  // prior
  theta_raw ~ normal(0, prior_sd);
  // likelihood, which we only evaluate conditionally
    y ~ categorical(softmax(theta));
generated quantities {
  vector[N] y_sim;
  for(i in 1:N) {
    y_sim[i] = categorical_rng(softmax(theta));

Now let’s compile and estimate this model with run_estimation=0, that is, in simulation mode. We’ll grab a random draw from the posterior, which will contain some fixed parameters, as well as the draws of fake data that we want.

compiled_model <- stan_model("categorical_model.stan")

sim_out <- sampling(compiled_model, data = list(N = 1000, 
                                                P = 5, 
                                                # Y should be real but fake data if we're simulating
                                                # new Y
                                                y = sample(1:2, 1000, replace = T), 
                                                run_estimation = 0,
                                                prior_sd = 100))


fake_data_matrix  <- sim_out %>% %>% 

summary_tbl <- apply(fake_data_matrix[1:5,], 1, summary)

Now we’ve got some fake data draws. Lets look at the distribution of outcomes for a very diffuse prior. Remember, each row of fake_data_matrix is associated with the corresponding prior draw.

1 2 3 4 5
Min. 2 2.000 3 2 4
1st Qu. 2 5.000 3 2 4
Median 2 5.000 3 2 4
Mean 2 4.985 3 2 4
3rd Qu. 2 5.000 3 2 4
Max. 2 5.000 3 2 4

The rows of the table are the relevant summary statistics of the simulated outcome for each draw. The columns correspond to the first five draws. As is clear, in each draw all the outcomes are precisely the same. This is because we’ve chosen a silly diffuse prior. (If this doesn’t make sense, check the definition of the softmax function and imagine what happens if you draw one \(\theta\) with a value of 400 or something). Let’s try it again with a more sensible prior.

sim_out_narrowprior <- sampling(compiled_model, data = list(N = 1000, 
                                                P = 5, 
                                                # Y should be real but fake data if we're simulating
                                                # new Y
                                                y = sample(1:2, 1000, replace = T), 
                                                run_estimation = 0,
                                                prior_sd = 1))

fake_data_matrix_2  <- sim_out_narrowprior %>% %>% 

summary_tbl_2 <- apply(fake_data_matrix_2[1:5,], 1, summary)
1 2 3 4 5
Min. 1.000 1.000 1.000 1.000 1.000
1st Qu. 2.000 3.000 2.000 2.000 3.000
Median 3.000 3.000 2.000 2.000 3.000
Mean 2.807 3.186 2.376 3.014 3.191
3rd Qu. 4.000 4.000 2.000 4.000 4.000
Max. 5.000 5.000 5.000 5.000 5.000

Much better! With a narrower prior we’re actually getting a variation in the outcome. Let’s grab our thetas and fake data from that run for a given draw. Say, 21. Then we’ll run the model again to make sure we can recapture those known thetas.

draw <- 21
y_sim <- extract(sim_out_narrowprior, pars = "y_sim")[[1]][draw,]
true_thetas <- extract(sim_out_narrowprior, pars = "theta_raw")[[1]][draw,]

# Now estimate the model with this "true" data and known parameters 

dgp_recapture_fit <- sampling(compiled_model, data = list(N = 1000, 
                                                P = 5, 
                                                # Note we now use y_sim
                                                y = y_sim, 
                                                # Note we not switch on run_estimation
                                                run_estimation = 1,
                                                prior_sd = 1))

How did we go?

Great! We just simulated some fake data from a Stan model, used it to discover that our original priors were insane, simulated again, and used the fake data to validate our model. Believe me—if you build models for a living, getting into the habit of doing this will save a lot of heartache.