Built-in way to plot prior and posterior distributions in RStan?

Is there a built-in way to plot prior and posterior distributions in the same plot? e.g. adding the prior for beta to the below?

stan_hist(fit, pars = c("beta"))

No, mostly because we do not know what the priors were for a general Stan program. We did it in the rstanarm::posterior_vs_prior function, which you could mostly copy for any particular Stan program where you have draws from the prior being used.

1 Like

Thanks so much ! Just to confirm I understand, for a general Stan program there is no direct way to draw from the prior (predictive distribution)?

You can do it within the same stan model, but it needs to be disconnected from the data.
Like, you can have the model defined as usual (e.g., with data y and parameter theta), then in generated quantities just generate random parameters (e.g., with y_prior and theta_prior), each according to their distributions.

E.g.,:

data {
int N;
vector[N] y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
mu ~ normal(0,2);
sigma ~ normal(0,2);
y ~ normal(mu, sigma);
}
generated quantities {
real mu_prior = normal_rng(0,2);
real sigma_prior = fabs(normal_rng(0,2)); // Or however you want to get a half-normal.
vector[N] y_prior;
for(n in 1:N){
y_prior[n] = normal_rng(mu_prior, sigma_prior);
}
}

I just wrote that on the fly; so there may be errors, and almost certainly better ways to do that. (I can’t recall if RNG is vectorized now, or whether you can directly sample a half-normal now). You could also just do all that in the parameters/model as well.

cool, thanks !! So it requires rewriting the distributions, e.g. normal(0,2) --> normal_rng(0,2) for mu?

If it is within generated quantities, yes.

Because when you see, e.g., y ~ normal(mu, sigma), that is the same thing as target += normal_lpdf(y | mu, sigma);, which is not relevant to generated quantities.

If you want to sample within the model itself, then you could do:

data {
...
}
parameters {
real mu;
real mu_prior;
real<lower=0> sigma;
real<lower=0> sigma_prior;
vector[N] y_prior;
}
model {
mu ~ normal(0,2);
mu_prior ~ normal(0,2);
sigma ~ normal(0,2);
sigma_prior ~ normal(0,2);
y ~ normal(mu, sigma);
y_prior ~ normal(mu_prior, sigma_prior);
}

But that adds complexity to the model. Conceptually, that is the same as generated quantities, but now Stan has a much higher dimensional space; it sees that as one joint model, although it is effectively two independent models (one with data, one without). That makes sampling harder overall, I believe, for no real gain over using generated quantities.