Calculate the values that greta arrays would take, given temporary values for the greta arrays on which they depend, and return them as numeric R arrays. This can be used to check the behaviour of your model or make predictions to new data after model fitting.

calculate(target, values = list(), precision = c("double", "single"))

Arguments

target

a greta array for which to calculate the value

values

a named list giving temporary values of the greta arrays with which target is connected, or an mcmc.list object returned by mcmc.

precision

the floating point precision to use when calculating values.

Value

A numeric R array with the same dimensions as target, giving the values it would take conditioned on the fixed values given by values.

Details

The greta arrays named in values need not be variables, they can also be other operations or even data.

At present, if values is a named list it must contain values for all of the variable greta arrays with which target is connected, even values are given for intermediate operations, or the target doesn't depend on the variable. That may be relaxed in a future release.

Examples

# NOT RUN {
# define a variable greta array, and another that is calculated from it
# then calculate what value y would take for different values of x
x <- normal(0, 1, dim = 3)
a <- lognormal(0, 1)
y <- sum(x ^ 2) + a
calculate(y, list(x = c(0.1, 0.2, 0.3), a = 2))

# if the greta array only depends on data,
# you can pass an empty list to values (this is the default)
x <- ones(3, 3)
y <- sum(x)
calculate(y)

# define a model
alpha <- normal(0, 1)
beta <- normal(0, 1)
sigma <- lognormal(1, 0.1)
mu <- alpha + iris$Petal.Length * beta
distribution(iris$Petal.Width) <- normal(mu, sigma)
m <- model(alpha, beta, sigma)

# calculate intermediate greta arrays, given some parameter values
calculate(mu[1:5], list(alpha = 1, beta = 2, sigma = 0.5))
calculate(mu[1:5], list(alpha = -1, beta = 0.2, sigma = 0.5))


# fit the model then calculate samples at a new greta array
draws <- mcmc(m, n_samples = 500)
petal_length_plot <- seq(min(iris$Petal.Length),
                         max(iris$Petal.Length),
                         length.out = 100)
mu_plot <- alpha + petal_length_plot * beta
mu_plot_draws <- calculate(mu_plot, draws)

# plot the draws
mu_est <- colMeans(mu_plot_draws[[1]])
plot(mu_est ~ petal_length_plot, type = "n",
     ylim = range(mu_plot_draws[[1]]))
apply(mu_plot_draws[[1]], 1, lines,
      x = petal_length_plot, col = grey(0.8))
lines(mu_est ~ petal_length_plot, lwd = 2)
# }