Pystan sampling from the posterior predictive

Interfaces PyStan
1

I love how clean and expressive the pystan api is.

sample = sm.sampling(data=data).to_dataframe()

I was just wondering if there’s a way to have the sampler sample from the likelihood conditional on each simulated parameter as well, in order to get posterior predictive observations added to this dataframe?

1 Like
2

You need to calculate them in generated quantities -block.

So for a simple model

data {
    int<lower=1> N;
    vector[N] y;
parameters {
    real mu;
model {
    y ~ normal(mu, 1);
}
generated quantities {
    vector[N] yhat;
    for (n in 1:N) {
        yhat[n] = normal_rng(mu, 1);
    }
}
1 Like
3

Perfect, thanks! Are vector expressions valid here too?

4

Some _rng functions are vectorised, so you don’t need an explicit loop.

1 Like
5

I have a question reading the line, not sure it is just an example code. I understand the posterior predictive distribution has the same mean as the posterior distribution, but it shouldn’t have the same variance.

Thanks!

6

It’s hardcoded to the model (so yeah, this is just an example model).

To fit it against the data, you need to define a parameter (sigma for example) and the use that (with the same structure as mu).

7

Here’s a slightly more complex model based on the original answer that does univariate OLS with prediction. In addition to the N observed datapoints you pass in for the regression fit, you also pass in K additional points that you’d like predictions for (they can be the same but don’t have to be). Then you can get point estimates via pystan.optimizing , or simulated observations via pystan.sampler that can be used to form Bayesian prediction intervals.

data {
    int<lower=0> N; // number of input
    vector[N] x;
    vector[N] y;
    int<lower=0> K; // number of predicted values
    vector[K] x_pred; // values we want predictions for
}
parameters {
    real alpha;
    real beta;
    real<lower=0> sigma;
}
model {
    y ~ normal(alpha +beta * x, sigma);
}
generated quantities {
    vector[K] yhat;
    for (n in 1:K) {
        yhat[n] = normal_rng(alpha +beta * x_pred[n], sigma);
    }
}
8

Hi, Sorry to jump on this so late.

Should yhat now also somehow contain many values for each to-be predicted value? I am only starting out with stan, but yhat[n] = normal_rng(alpha +beta * x_pred[n], sigma); looks like there is a single value predicted.

9

Yes. Only single value is predicted here (per draw). You could add another dimension and predict more values (per draw)