Generating Quantities of Interest from a Fitted Model

The generate_quantities method allows you to generate additional quantities of interest from a fitted model without re-running the sampler. Instead, you write a modified version of the original Stan program and add a generated quantities block or modify the existing one which specifies how to compute the new quantities of interest. Running the generate_quantities method on the new program together with sampler outputs (i.e., a set of draws) from the fitted model runs the generated quantities block of the new program using the the existing sample by plugging in the per-draw parameter estimates for the computations in the generated quantities block.

This method requires sub-argument fitted_params which takes as its value an existing Stan CSV file that contains a parameter values from an equivalent model, i.e., a model with the same parameters block, conditioned on the same data.

The generated quantities block computes quantities of interest (QOIs) based on the data, transformed data, parameters, and transformed parameters. It can be used to:

  • generate simulated data for model testing by forward sampling
  • generate predictions for new data
  • calculate posterior event probabilities, including multiple comparisons, sign tests, etc.
  • calculate posterior expectations
  • transform parameters for reporting
  • apply full Bayesian decision theory
  • calculate log likelihoods, deviances, etc. for model comparison

For an overview of the uses of this feature, see the Stan User’s Guide section on Stand-alone generated quantities and ongoing prediction.

Example

To illustrate how this works we use the generate_quantities method to do posterior predictive checks using the estimate of theta given the example bernoulli model and data, following the posterior predictive simulation procedure in the Stan User’s Guide.

We write a program bernoulli_ppc.stan which contains the following generated quantities block, with comments to explain the procedure:

generated quantities {
  array[N] int y_sim;
  // use current estimate of theta to generate new sample
  for (n in 1:N) {
    y_sim[n] = bernoulli_rng(theta);
  }
  // estimate theta_rep from new sample
  real<lower=0, upper=1> theta_rep = sum(y_sim) * 1.0 / N;
}

The rest of the program is the same as in bernoulli.stan.

The generate_method requires the sub-argument fitted_params which takes as its value the name of a Stan CSV file. The per-draw parameter values from the fitted_params file will be used to run the generated quantities block.

If we run the bernoulli.stan program for a single chain to generate a sample in file bernoulli_fit.csv:

> ./bernoulli sample data file=bernoulli.data.json output file=bernoulli_fit.csv

Then we can run the bernoulli_ppc.stan to carry out the posterior predictive checks:

> ./bernoulli_ppc generate_quantities fitted_params=bernoulli_fit.csv \
                  data file=bernoulli.data.json \
                  output file=bernoulli_ppc.csv

The output file bernoulli_ppc.csv contains only the values for the variables declared in the generated quantities block, i.e., theta_rep and the elements of y_sim:

# model = bernoulli_ppc_model
# method = generate_quantities
#   generate_quantities
#     fitted_params = bernoulli_fit.csv
# id = 1 (Default)
# data
#   file = bernoulli.data.json
# init = 2 (Default)
# random
#   seed = 2983956445 (Default)
# output
#   file = output.csv (Default)
y_sim.1,y_sim.2,y_sim.3,y_sim.4,y_sim.5,y_sim.6,y_sim.7,y_sim.8,y_sim.9,y_sim.10,theta_rep
1,1,1,0,0,0,1,1,0,1,0.6
1,1,0,1,0,0,1,0,1,0,0.5
1,0,1,1,1,1,1,1,0,1,0.8
0,1,0,1,0,1,0,1,0,0,0.4
1,0,0,0,0,0,0,0,0,0,0.1
0,0,0,0,0,1,1,1,0,0,0.3
0,0,1,0,1,0,0,0,0,0,0.2
1,0,1,0,1,1,0,1,1,0,0.6
...

Given the current implementation, to see the fitted parameter values for each draw, create a copy variable in the generated quantities block, e.g.:

generated quantities {
  array[N] int y_sim;
  // use current estimate of theta to generate new sample
  for (n in 1:N) {
    y_sim[n] = bernoulli_rng(theta);
  }
  real<lower=0, upper=1> theta_cp = theta;
  // estimate theta_rep from new sample
  real<lower=0, upper=1> theta_rep = sum(y_sim) * 1.0 / N;
}

Now the output is slightly more interpretable: theta_cp is the same as the theta used to generate the values y_sim[1] through y_sim[1]. Comparing columns theta_cp and theta_rep allows us to see how the uncertainty in our estimate of theta is carried forward into our predictions:

y_sim.1,y_sim.2,y_sim.3,y_sim.4,y_sim.5,y_sim.6,y_sim.7,y_sim.8,y_sim.9,y_sim.10,theta_cp,theta_rep
0,1,1,0,1,0,0,1,1,0,0.545679,0.5
1,1,1,1,1,1,0,1,1,0,0.527164,0.8
1,1,1,1,0,1,1,1,1,0,0.529116,0.8
1,0,1,1,1,1,0,0,1,0,0.478844,0.6
0,1,0,0,0,0,1,0,1,0,0.238793,0.3
0,0,0,0,0,1,1,0,0,0,0.258294,0.2
1,1,1,0,0,0,0,0,0,0,0.258465,0.3

Errors

The fitted_params file must be a Stan CSV file; attempts to use a regular CSV file will result an error message of the form:

Error reading fitted param names from sample csv file <filename.csv>

The fitted_params file must contain columns corresponding to legal values for all parameters defined in the model. If any parameters are missing, the program will exit with an error message of the form:

Error reading fitted param names from sample csv file <filename.csv>

The parameter values of the fitted_params are on the constrained scale and must obey all constraints. For example, if we modify the contents of the first reported draw in bernoulli_fit.csv so that the value of theta is outside the declared bounds real<lower=0, upper=1>, the program will return the following error message:

Exception: lub_free: Bounded variable is 1.21397, but must be in the interval [0, 1] \
(in 'bernoulli_ppc.stan', line 5, column 2 to column 30)
Back to top