Two more things:
1a. We can do extract(fit)$theta to get that n_sims x 2 x 3 x 4 array. I’d also thought it could make sense to do it as extract(fit, “theta”). But for some reason now when we do it that way in R, we have to do extract(fit, “theta”)$theta. That is, extract returns a list, and even fi only one named parameter is extracted, it still makes it a list with one item; hence as a user I need to do extract(fit, “theta”)$theta or extract(fit, “theta”)[[1]]. Maybe that’s just the way R is, it’s no big deal.
1b. I’d also like extract_one_draw(fit, i) which returns the list of all the parameters for one draw (the i’th draw). For example, if alpha is a scalar, beta is a vector of length 10, and theta is a 2 x 3 x 4 array, then
extract_one_draw(fit, 34)$alpha is a scalar
extract_one_draw(fit, 34)$beta is a vector of length 10
extract_one_draw(fit, 34)$theta is a 2 x 3 x 4 array,
all corresponding to the same simulation draw.
These are equivalent to extract(fit)$alpha[34], extract(fti)$beta[34,], and extract(fit)$theta[34,], but the point is that we can grab them all at once. This is super-important for propagating uncertainty and for probabilistic programming in general.