Thanks! Yeah, that would probably be the most common use case. Though I worry that it would be non-intuitive given other functions operating on models (mcmc
, opt
, etc.) always return summaries of the greta arrays named in the model call. Not sure which is best.
I’ve been thinking again about the marginalisation interface in the feature/marginalise
branch, and also doing predictive simulations (when you get new data to predict to after model fitting). In both of those use cases, you might want to create a new variable greta array after creating the model and doing MCMC, and then simulate their values given the MCMC samples (like you can already do for deterministic greta arrays with calculate()
). You can’t do that with the current interface, since you have to pass in the model object, which wouldn’t have the new greta arrays.
In light of that, maybe having greta arrays be the simulation targets would be better, like in @cboettig’s post. so then the interface would look like this (no model passed in):
simulate(targets, nsim = 1, seed = NULL, values = list(), precision = c("double", "single"), ...)
and you’d use it like this:
# build a greta model
n <- 10
y <- rnorm(n)
y <- as_data(y)
library(greta)
sd <- lognormal(1, 2)
mu <- normal(0, 1, dim = n)
distribution(y) <- normal(mu, sd)
m <- model(mu, sd)
# prior simulation:
# simulate one random draw of mu and sd from the model prior:
sims <- simulate(list(mu, sd))
# 100 simulations of mu and sd)
sims <- simulate(list(mu, sd), nsim = 100)
# 100 simulations of y and mu
# NB: y must be a greta array here
sims <- simulate(list(y, mu), nsim = 100)
# conditional prior simulation:
# 100 simulations of y and mu, with sd fixed at 3
sims <- simulate(list(y, mu), nsim = 100, values = list(sd = 3))
# posterior simulation:
# generate posterior samples from the model as usual
draws <- mcmc(m)
# draw 100 simulations of y, conditioned on 100 random posterior
# samples of all the model parameters
sims <- simulate(y, nsim = 100, values = draws)
This way you would have to explicitly pass in the things you want to simulate. It wouldn’t be able to guess what you want with defaults, but at least you wouldn’t be confused!
A downside of this is that it doesn’t match the pattern for the simulate
S3 class, where the first argument is a (fitted) model object.
So maybe it should be renamed to something else (sample
is already taken) or maybe this functionality should be rolled into the same function as calculate
(enabling a list of targets in that function)?
Goddamnit, sane API design is much harder than writing code!