Just to be clear on the specs:

I will use the term “parameter” to refer to a named parameter. For example, if “a” is a vector of length 10, then “a” is one parameter. I would not refer to “a[3]” as a parameter but rather as an element of the parameter “a”. I think the integrity of the named parameter is a key part of Stan being a probabilistic programming language.

Suppose we fit a Stan model with 3 parameters:

parameters {

vector[10] a;

real b[3,5];

real<lower=0> sigma;

}

We fit the model and get a stanfit object, which I’ll call “fit”.

So here’s what I’d like:

We’d have a function called “extract” (or something like that) which would take 2 arguments: the stanfit object and an instruction of what to pull out. Thus, extract <- function(stanfit, operation, …), where … are arguments to operation.

For example:

foo <- extract(fit, operation=sample, 44)

would return a list of length 3. The first element of the list would be named “a” and would be a vector of length 10, with the values of the 44th draw of “a” from the posterior. The second element of the list would be named “b” and would be an array of dimension [3,5], etc.

For another example,

foo <- extract(fit, operation=“median”)

would return a list of length 3. The first element of the list would be named “a” and would be a vector of length 10, corresponding to the elementwise posterior median of this parameter. The second element of the list would be named “b” and would be an array of dimension [3,5], etc.

Or:

foo <- extract(fit, operation=“quantile”, 0.80)

would return a list where each item is the elementwise 80% quantile from the posterior simulation.

I’m no expert on what is good R practice. Perhaps it would make more sense to have two functions, one for posterior draws and another for posterior summaries. It could go either way.

One thing that seems clear, though, is that everything would be done in an elementwise way. This is the “probabilistic programming” thing where, after fitting the Stan model, we can directly grab items that are the same size and shape as the “a”, “b”, and “sigma” in the model.

This has been coming up in just about every example I’ve been doing.

It looks like from Ben’s last reply that there might be a way to hack together this function using relists of some sort? If so, I think it would be a good idea for us to create the function as then we can all use it while we await cleaner future versions.