Prior Predictive Checks with RStan

Hi everybody,

I would like to perform a prior predictive check for my Stan model in RStan, following https://mc-stan.org/docs/2_24/stan-users-guide/prior-predictive-checks.html

My Stan Code:

data {
  int<lower=0> N;
  vector[N] x_1;
}
generated quantities {
  real<lower=0> alpha = normal_rng(0, 0.5);
  real beta_1 = normal_rng(0.1, 1);
  real<lower=0> sigma;
  real y_sim[N] = normal_rng(alpha + beta_1*x_1, sigma);
}

My R Code:

x_1 = sample(c(0,40), 1000, replace = TRUE)  # reflecting true data 
model_data = list(N = length(x_1),
                   x_1)

fit = stan(file = "scripts/simple_regression_prior_check.stan",
           data = model_data,
           chains = 4,
           warmup = 1000,
           iter = 2000,
           cores =4)

I get the following error:
*Error in mod$fit_ptr() : *

  • Exception: variable does not exist; processing stage=data initialization; variable name=x_1; base type=vector_d (in ‘modelafb3458c0570_simple_regression_prior_check’ at line 3)*

failed to create the sampler; sampling not done

I am fairly new to stan, so apologies if it is obvious that this is not the way to do it, but I have not found any ressources on how to properly do a prior posterior check in rstan.

I know that it is seemingly also possible to do this with brms or bayesplot, but I have not yet advanced to these package and as the model is fairly simple, I would like to do it without or at least understand why what I am doing does not work.

Any thoughts and comments appreciated.

1 Like

Each element in the data list should be named. Try the following:

model_data <- list(N=length(x_1), x_1=x_1)
1 Like

Yes, thank you. I did that, now I get

Stan model 'simple_regression_prior_check' does not contain samples

I am in general asking if this is the way to calculate prior predictive checks in RStan?

It appears there is additional problem with your model. Try running with chains = 1, this should show all messages directly in your console.

I can see two problems with your code:

real<lower=0> alpha = normal_rng(0, 0.5);

normal_rng will happily return negative values, violating the constraint you gave (slightly counter-intuitively, constraints behave differently in parameters block and elsewhere - outside of paramters they are simply something the program checks for and doesn’t enforce them, so you would need fabs(normal_rng(0, 0.5)) to satisfy the constraint.

sigma is not initialized before it is used…

Best of luck with your checks!

1 Like