I’d like to have fit$theta because this comes up a lot! I’d also like to do something like foo(fit$theta, median) or foo(fit, “theta”, median) that would grab the median over the simulation draws (or it could be foo(fit, “theta”, mean) or foo(fit, “theta”, quantile, c(0.25, 0.75)) etc. Currently I waste a lot of time in R grabbing simulations, then trying to figure out how to use “apply” or whatever to get summaries. This is a particular issue if theta is a vector or matrix or array. What I’d generally like is pointwise summaries; e.g., if theta is a 2 x 3 x 4 array, then fit$theta would return a n_sims x 2 x 3 x 4 array, and foo(fit, “theta”, median) would return a 2 x 3 x 4 array.
I’m also thinking the default would be to concatenate all chains together, but of course there should also be an option to grab one chain at a time, which could come up for diagnostics or whatever.