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

Interfaces RStan
1

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"))

Rplot
Rplot.png683×549 6.09 KB

2

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
3

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)?

4

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.

1 Like
5

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

6

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.

1 Like