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"))



a greta array for which to calculate the value


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


the floating point precision to use when calculating values.


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


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.


# 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)

# 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),
                         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)
# }