I would suggest asking @andrewgelman, who often has opinions on interface issues.
On the specific topic of fit$extract(...)
as extract(fit, ...)
, they’re interchangeable to me.
On the other hand, fit$theta
feels odd to me as it feels like the structure of the fit
object changes.
I’m all trial and error and str()
in R, so it doesn’t matter that much to me.
As for StanDimensions x Chains x MainIterations
, yes, that makes most sense to me (assuming we mean Cstyle indexing here, where we first provide a Stan dimension (picking out a single scalar parameter) then get a collection of chains for that parameter). It lets you get the chain x iterations for each parameter dimension, which is what we need for the posterior.
I almost never look at my own chains, so it’d be nice to also have the collapsed option of getting StanDimensions x Draws
. It doesn’t need to be permuted.
Is the Stan fit object in R going to be like the current one in that it includes input, model binary, etc.?
What’s the purpose of flattening theta[2]
(a size 2 array or 2D vector)? I won’t ever have a use for that. What I think we need is be able to give a parameter (say in string form with indexes as it’s done now, say theta[1]
) and get either

N chains, each consisting of M draws of theta[1]
, or

N * M draws of theta[1]
.
I’ve never needed (1).
What I do find clumsy is having to use paste to manipulate strings:
print(fit, pars = c("alpha", paste("theta[", 1:5, "]"))