Package 'greta'

Title: Simple and Scalable Statistical Modelling in R
Description: Write statistical models in R and fit them by MCMC and optimisation on CPUs and GPUs, using Google 'TensorFlow'. greta lets you write your own model like in BUGS, JAGS and Stan, except that you write models right in R, it scales well to massive datasets, and it’s easy to extend and build on. See the website for more information, including tutorials, examples, package documentation, and the greta forum.
Authors: Nick Golding [aut] , Nicholas Tierney [aut, cre] , Simon Dirmeier [ctb], Adam Fleischhacker [ctb], Shirin Glander [ctb], Martin Ingram [ctb], Lee Hazel [ctb], Lionel Hertzog [ctb], Tiphaine Martin [ctb], Matt Mulvahill [ctb], Michael Quinn [ctb], David Smith [ctb], Paul Teetor [ctb], Jian Yen [ctb]
Maintainer: Nicholas Tierney <[email protected]>
License: Apache License 2.0
Version: 0.4.5
Built: 2024-07-09 05:05:32 UTC

Help Index

convert other objects to greta arrays


define an object in an R session as a data greta array for use as data in a greta model.





an R object that can be coerced to a greta_array (see details).


as_data() can currently convert R objects to greta_arrays if they are numeric or logical vectors, matrices or arrays; or if they are dataframes with only numeric (including integer) or logical elements. Logical elements are always converted to numerics. R objects cannot be converted if they contain missing (NA) or infinite (-Inf or Inf) values.


## Not run: 

# numeric/integer/logical vectors, matrices and arrays can all be coerced to
# data greta arrays

vec <- rnorm(10)
mat <- matrix(seq_len(3 * 4), nrow = 3)
arr <- array(sample(c(TRUE, FALSE), 2 * 2 * 2, replace = TRUE),
  dim = c(2, 2, 2)
(a <- as_data(vec))
(b <- as_data(mat))
(c <- as_data(arr))

# dataframes can also be coerced, provided all the columns are numeric,
# integer or logical
df <- data.frame(
  x1 = rnorm(10),
  x2 = sample(1L:10L),
  x3 = sample(c(TRUE, FALSE), 10, replace = TRUE)
(d <- as_data(df))

## End(Not run)

calculate greta arrays given fixed values


Calculate the values that greta arrays would take, given temporary, or simulated values for the greta arrays on which they depend. This can be used to check the behaviour of your model, make predictions to new data after model fitting, or simulate datasets from either the prior or posterior of your model.


  values = list(),
  nsim = NULL,
  seed = NULL,
  precision = c("double", "single"),
  trace_batch_size = 100



one or more greta_arrays for which to calculate the value


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


an optional positive integer scalar for the number of responses to simulate if stochastic greta arrays are present in the model - see Details.


an optional seed to be used in set.seed immediately before the simulation so as to generate a reproducible sample


the floating point precision to use when calculating values.


the number of posterior samples to process at a time when target is a greta_mcmc_list object; reduce this to reduce memory demands


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.

If the model contains stochastic greta arrays; those with a distribution, calculate can be used to sample from these distributions (and all greta arrays that depend on them) by setting the nsim argument to a positive integer for the required number of samples. If values is specified (either as a list of fixed values or as draws), those values will be used, and remaining variables will be sampled conditional on them. Observed data with distributions (i.e. response variables defined with distribution() can also be sampled, provided they are defined as greta arrays. This behaviour can be used for a number of tasks, like simulating datasets for known parameter sets, simulating parameters and data from a set of priors, or simulating datasets from a model posterior. See some examples of these below.


Values of the target greta array(s), given values of the greta arrays on which they depend (either specified in values or sampled from their priors). If values is a greta_mcmc_list() and nsim = NULL, this will be a greta_mcmc_list object of posterior samples for the target greta arrays. Otherwise, the result will be a named list of numeric R arrays. If nsim = NULL the dimensions of returned numeric R arrays will be the same as the corresponding greta arrays, otherwise an additional dimension with nsim elements will be prepended, to represent multiple simulations.


## 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, values = list(x = c(0.1, 0.2, 0.3), a = 2))

# by setting nsim, you can also sample values from their priors
calculate(y, nsim = 3)

# you can combine sampling and fixed values
calculate(y, values = list(a = 2), nsim = 3)

# 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)
y <- as_data(iris$Petal.Width)
mu <- alpha + iris$Petal.Length * beta
distribution(y) <- normal(mu, sigma)
m <- model(alpha, beta, sigma)

# sample values of the parameters, or different observation data (y), from
# the priors (useful for prior # predictive checking) - see also
# ?simulate.greta_model
calculate(alpha, beta, sigma, nsim = 100)
calculate(y, nsim = 100)

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

# simulate datasets given fixed parameter values
calculate(y, values = list(alpha = -1, beta = 0.2, sigma = 0.5), nsim = 10)

# you can use calculate in conjunction with posterior samples from MCMC, e.g.
# sampling different observation datasets, given a random set of these
# posterior samples - useful for posterior predictive model checks
draws <- mcmc(m, n_samples = 500)
calculate(y, values = draws, nsim = 100)

# you can use calculate on greta arrays created even after the inference on
# the model - e.g. to plot response curves
petal_length_plot <- seq(min(iris$Petal.Length),
  length.out = 100
mu_plot <- alpha + petal_length_plot * beta
mu_plot_draws <- calculate(mu_plot, values = 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)

# trace_batch_size can be changed to trade off speed against memory usage
# when calculating. These all produce the same result, but have increasing
# memory requirements:
mu_plot_draws_1 <- calculate(mu_plot,
  values = draws,
  trace_batch_size = 1
mu_plot_draws_10 <- calculate(mu_plot,
  values = draws,
  trace_batch_size = 10
mu_plot_draws_inf <- calculate(mu_plot,
  values = draws,
  trace_batch_size = Inf

## End(Not run)

Cholesky Factor to Symmetric Matrix


Evaluate ⁠t(x) \%*\% x⁠ efficiently, where x is the (upper-triangular) Cholesky factor of a symmetric, positive definite square matrix. I.e. it is the inverse of chol.





a square, upper triangular matrix representing the Cholesky factor of a symmetric, positive definite square matrix


# a symmetric, positive definite square matrix
y <- rWishart(1, 4, diag(3))[, , 1]
u <- chol(y)
identical(y, chol2symm(u))
identical(chol2symm(u), t(u) %*% u)
## Not run: 
u_greta <- cholesky_variable(3)
y_greta <- chol2symm(u)

## End(Not run)

define a distribution over data


distribution defines probability distributions over observed data, e.g. to set a model likelihood.


distribution(greta_array) <- value




a data greta array. For the assignment method it must not already have a probability distribution assigned


a greta array with a distribution (see distributions())


The extract method returns the greta array if it has a distribution, or NULL if it doesn't. It has no real use-case, but is included for completeness


## Not run: 

# define a model likelihood

# observed data and mean parameter to be estimated
# (explicitly coerce data to a greta array so we can refer to it later)
y <- as_data(rnorm(5, 0, 3))

mu <- uniform(-3, 3)

# define the distribution over y (the model likelihood)
distribution(y) <- normal(mu, 1)

# get the distribution over y

## End(Not run)

probability distributions


These functions can be used to define random variables in a greta model. They return a variable greta array that follows the specified distribution. This variable greta array can be used to represent a parameter with prior distribution, combined into a mixture distribution using mixture(), or used with distribution() to define a distribution over a data greta array.


uniform(min, max, dim = NULL)

normal(mean, sd, dim = NULL, truncation = c(-Inf, Inf))

lognormal(meanlog, sdlog, dim = NULL, truncation = c(0, Inf))

bernoulli(prob, dim = NULL)

binomial(size, prob, dim = NULL)

beta_binomial(size, alpha, beta, dim = NULL)

negative_binomial(size, prob, dim = NULL)

hypergeometric(m, n, k, dim = NULL)

poisson(lambda, dim = NULL)

gamma(shape, rate, dim = NULL, truncation = c(0, Inf))

inverse_gamma(alpha, beta, dim = NULL, truncation = c(0, Inf))

weibull(shape, scale, dim = NULL, truncation = c(0, Inf))

exponential(rate, dim = NULL, truncation = c(0, Inf))

pareto(a, b, dim = NULL, truncation = c(0, Inf))

student(df, mu, sigma, dim = NULL, truncation = c(-Inf, Inf))

laplace(mu, sigma, dim = NULL, truncation = c(-Inf, Inf))

beta(shape1, shape2, dim = NULL, truncation = c(0, 1))

cauchy(location, scale, dim = NULL, truncation = c(-Inf, Inf))

chi_squared(df, dim = NULL, truncation = c(0, Inf))

logistic(location, scale, dim = NULL, truncation = c(-Inf, Inf))

f(df1, df2, dim = NULL, truncation = c(0, Inf))

multivariate_normal(mean, Sigma, n_realisations = NULL, dimension = NULL)

wishart(df, Sigma)

lkj_correlation(eta, dimension = 2)

multinomial(size, prob, n_realisations = NULL, dimension = NULL)

categorical(prob, n_realisations = NULL, dimension = NULL)

dirichlet(alpha, n_realisations = NULL, dimension = NULL)

dirichlet_multinomial(size, alpha, n_realisations = NULL, dimension = NULL)


min, max

scalar values giving optional limits to uniform variables. Like lower and upper, these must be specified as numerics, they cannot be greta arrays (though see details for a workaround). Unlike lower and upper, they must be finite. min must always be less than max.


the dimensions of the greta array to be returned, either a scalar or a vector of positive integers. See details.

mean, meanlog, location, mu

unconstrained parameters

sd, sdlog, sigma, lambda, shape, rate, df, scale, shape1, shape2, alpha, beta, df1, df2, a, b, eta

positive parameters, alpha must be a vector for dirichlet and dirichlet_multinomial.


a length-two vector giving values between which to truncate the distribution, similarly to the lower and upper arguments to variable()


probability parameter (⁠0 < prob < 1⁠), must be a vector for multinomial and categorical

size, m, n, k

positive integer parameter


positive definite variance-covariance matrix parameter


the number of independent realisation of a multivariate distribution


the dimension of a multivariate distribution


The discrete probability distributions (bernoulli, binomial, negative_binomial, poisson, multinomial, categorical, dirichlet_multinomial) can be used when they have fixed values (e.g. defined as a likelihood using distribution(), but not as unknown variables.

For univariate distributions dim gives the dimensions of the greta array to create. Each element of the greta array will be (independently) distributed according to the distribution. dim can also be left at its default of NULL, in which case the dimension will be detected from the dimensions of the parameters (provided they are compatible with one another).

For multivariate distributions (multivariate_normal(), multinomial(), categorical(), dirichlet(), and dirichlet_multinomial()) each row of the output and parameters corresponds to an independent realisation. If a single realisation or parameter value is specified, it must therefore be a row vector (see example). n_realisations gives the number of rows/realisations, and dimension gives the dimension of the distribution. I.e. a bivariate normal distribution would be produced with multivariate_normal(..., dimension = 2). The dimension can usually be detected from the parameters.

multinomial() does not check that observed values sum to size, and categorical() does not check that only one of the observed entries is 1. It's the user's responsibility to check their data matches the distribution!

The parameters of uniform must be fixed, not greta arrays. This ensures these values can always be transformed to a continuous scale to run the samplers efficiently. However, a hierarchical uniform parameter can always be created by defining a uniform variable constrained between 0 and 1, and then transforming it to the required scale. See below for an example.

Wherever possible, the parameterisations and argument names of greta distributions match commonly used R functions for distributions, such as those in the stats or extraDistr packages. The following table states the distribution function to which greta's implementation corresponds:

greta reference
uniform stats::dunif
normal stats::dnorm
lognormal stats::dlnorm
bernoulli extraDistr::dbern
binomial stats::dbinom
beta_binomial extraDistr::dbbinom
negative_binomial stats::dnbinom
hypergeometric stats::dhyper
poisson stats::dpois
gamma stats::dgamma
inverse_gamma extraDistr::dinvgamma
weibull stats::dweibull
exponential stats::dexp
pareto extraDistr::dpareto
student extraDistr::dlst
laplace extraDistr::dlaplace
beta stats::dbeta
cauchy stats::dcauchy
chi_squared stats::dchisq
logistic stats::dlogis
f stats::df
multivariate_normal mvtnorm::dmvnorm
multinomial stats::dmultinom
categorical stats::dmultinom (size = 1)
dirichlet extraDistr::ddirichlet
dirichlet_multinomial extraDistr::ddirmnom
wishart stats::rWishart
lkj_correlation rethinking::dlkjcorr


## Not run: 

# a uniform parameter constrained to be between 0 and 1
phi <- uniform(min = 0, max = 1)

# a length-three variable, with each element following a standard normal
# distribution
alpha <- normal(0, 1, dim = 3)

# a length-three variable of lognormals
sigma <- lognormal(0, 3, dim = 3)

# a hierarchical uniform, constrained between alpha and alpha + sigma,
eta <- alpha + uniform(0, 1, dim = 3) * sigma

# a hierarchical distribution
mu <- normal(0, 1)
sigma <- lognormal(0, 1)
theta <- normal(mu, sigma)

# a vector of 3 variables drawn from the same hierarchical distribution
thetas <- normal(mu, sigma, dim = 3)

# a matrix of 12 variables drawn from the same hierarchical distribution
thetas <- normal(mu, sigma, dim = c(3, 4))

# a multivariate normal variable, with correlation between two elements
# note that the parameter must be a row vector
Sig <- diag(4)
Sig[3, 4] <- Sig[4, 3] <- 0.6
theta <- multivariate_normal(t(rep(mu, 4)), Sig)

# 10 independent replicates of that
theta <- multivariate_normal(t(rep(mu, 4)), Sig, n_realisations = 10)

# 10 multivariate normal replicates, each with a different mean vector,
# but the same covariance matrix
means <- matrix(rnorm(40), 10, 4)
theta <- multivariate_normal(means, Sig, n_realisations = 10)

# a Wishart variable with the same covariance parameter
theta <- wishart(df = 5, Sigma = Sig)

## End(Not run)

extract, replace and combine greta arrays


Generic methods to extract and replace elements of greta arrays, or to combine greta arrays.



a greta array

i, j

indices specifying elements to extract or replace


a single integer, as in utils::head() and utils::tail()

nrow, ncol

optional dimensions for the resulting greta array when x is not a matrix.


for ⁠[<-⁠ a greta array to replace elements, for ⁠dim<-⁠ either NULL or a numeric vector of dimensions


either further indices specifying elements to extract or replace ([), or multiple greta arrays to combine (cbind(), rbind() & c()), or additional arguments (rep(), head(), tail())

drop, recursive

generic arguments that are ignored for greta arrays


diag() can be used to extract or replace the diagonal part of a square and two-dimensional greta array, but it cannot be used to create a matrix-like greta array from a scalar or vector-like greta array. A static diagonal matrix can always be created with e.g. diag(3), and then converted into a greta array.

Also note that since R 4.0.0, head and tail methods for arrays changed to print a vector rather than maintain the array structure. The greta package supports both methods, and will do so based on which version of R you are using.


# extract
x[i, j, ..., drop = FALSE]
head(x, n = 6L, ...)
tail(x, n = 6L, ...)
diag(x, nrow, ncol)

# replace
x[i] <- value
x[i, j, ...] <- value
diag(x) <- value

# combine
c(..., recursive = FALSE)
rep(x, times, ..., recursive = FALSE)

# get and set dimensions
dim(x) <- value


## Not run: 

x <- as_data(matrix(1:12, 3, 4))

# extract and replace
x[1:3, ]
x[, 2:4] <- 1:9
e <- diag(x)
diag(x) <- e + 1

# combine
cbind(x[, 2], x[, 1])
rbind(x[1, ], x[3, ])
abind(x[1, ], x[3, ], along = 1)
c(x[, 1], x)
rep(x[, 2], times = 3)

## End(Not run)

functions for greta arrays


This is a list of functions (mostly from base R) that are currently implemented to transform greta arrays. Also see operators and transforms.


TensorFlow only enables rounding to integers, so round() will error if digits is set to anything other than 0.

Any additional arguments to chol(), chol2inv, and solve() will be ignored, see the TensorFlow documentation for details of these routines.

sweep() only works on two-dimensional greta arrays (so MARGIN can only be either 1 or 2), and only for subtraction, addition, division and multiplication.

tapply() works on column vectors (2D greta arrays with one column), and INDEX cannot be a greta array. Currently five functions are available, and arguments passed to ... are ignored.

cospi(), sinpi(), and tanpi() do not use the computationally more stable routines to compute cos(x * pi) etc. that are available in R under some operating systems. Similarly trigamma() uses TensorFlow's polygamma function, resulting in lower precision than R's equivalent.


 # logarithms and exponentials

 # miscellaneous mathematics

 # rounding of numbers
 round(x, digits = 0)

 # trigonometry

 # special mathematical functions
 choose(n, k)
 lchoose(n, k)

 # matrix operations
 chol(x, ...)
 chol2inv(x, ...)
 solve(a, b, ...)
 kronecker(X, Y, FUN = c('*', '/', '+', '-'))

 # reducing operations
 sum(..., na.rm = TRUE)
 prod(..., na.rm = TRUE)
 min(..., na.rm = TRUE)
 max(..., na.rm = TRUE)

 # cumulative operations

 # solve an upper or lower triangular system
 backsolve(r, x, k = ncol(r), upper.tri = TRUE,
           transpose = FALSE)
 forwardsolve(l, x, k = ncol(l), upper.tri = FALSE,
              transpose = FALSE)

 # miscellaneous operations
 aperm(x, perm)
 apply(x, MARGIN, FUN = c("sum", "max", "mean", "min",
                          "prod", "cumsum", "cumprod"))
 sweep(x, MARGIN, STATS, FUN = c('-', '+', '/', '*'))
 tapply(X, INDEX, FUN = c("sum", "max", "mean", "min", "prod"), ...)


## Not run: 

x <- as_data(matrix(1:9, nrow = 3, ncol = 3))
a <- log(exp(x))
b <- log1p(expm1(x))
c <- sign(x - 5)
d <- abs(x - 5)

z <- t(a)

y <- sweep(x, 1, e, "-")

## End(Not run)

greta: simple and scalable statistical modelling in R


greta lets you write statistical models interactively in native R code, then sample from them efficiently using Hamiltonian Monte Carlo.

The computational heavy lifting is done by TensorFlow, Google's automatic differentiation library. So greta is particularly fast where the model contains lots of linear algebra, and greta models can be run across CPU clusters or on GPUs.

See the simple example below, and take a look at the greta website for more information including tutorials and examples.


Maintainer: Nicholas Tierney [email protected] (ORCID)


Other contributors:

  • Simon Dirmeier [contributor]

  • Adam Fleischhacker [contributor]

  • Shirin Glander [contributor]

  • Martin Ingram [contributor]

  • Lee Hazel [contributor]

  • Lionel Hertzog [contributor]

  • Tiphaine Martin [contributor]

  • Matt Mulvahill [contributor]

  • Michael Quinn [contributor]

  • David Smith [contributor]

  • Paul Teetor [contributor]

  • Jian Yen [contributor]

See Also

Useful links:


## Not run: 
# a simple Bayesian regression model for the iris data

# priors
int <- normal(0, 5)
coef <- normal(0, 3)
sd <- lognormal(0, 3)

# likelihood
mean <- int + coef * iris$Petal.Length
distribution(iris$Sepal.Length) <- normal(mean, sd)

# build and sample
m <- model(int, coef, sd)
draws <- mcmc(m, n_samples = 100)

## End(Not run)

Retrieve python installation or error details.


These functions retrieve installation or error information output by python when running install_miniconda(), conda_create(), conda_install(), or when encountering a TensorFlow numerical problem.










## Not run: 

## End(Not run)

Greta Situation Report


This checks if Python, Tensorflow, Tensorflow Probability, and the greta conda environment are available, and also loads and initialises python




Message if greta is ready to use


## Not run: 

## End(Not run)

Statistical inference on greta models.


Carry out statistical inference on greta models by MCMC or likelihood/posterior optimisation.


  sampler = hmc(),
  n_samples = 1000,
  thin = 1,
  warmup = 1000,
  chains = 4,
  n_cores = NULL,
  verbose = TRUE,
  pb_update = 50,
  one_by_one = FALSE,
  initial_values = initials(),
  trace_batch_size = 100


  n_samples = 1000,
  thin = 1,
  n_cores = NULL,
  verbose = TRUE,
  pb_update = 50,
  one_by_one = FALSE,
  trace_batch_size = 100


  optimiser = bfgs(),
  max_iterations = 100,
  tolerance = 1e-06,
  initial_values = initials(),
  adjust = TRUE,
  hessian = FALSE



greta_model object


sampler used to draw values in MCMC. See samplers() for options.


number of MCMC samples to draw per chain (after any warm-up, but before thinning)


MCMC thinning rate; every thin samples is retained, the rest are discarded


number of samples to spend warming up the mcmc sampler (moving chains toward the highest density area and tuning sampler hyperparameters).


number of MCMC chains to run


the maximum number of CPU cores used by each sampler (see details).


whether to print progress information to the console


how regularly to update the progress bar (in iterations). If pb_update is less than or equal to thin, it will be set to thin + 1 to ensure at least one saved iteration per pb_update iterations.


whether to run TensorFlow MCMC code one iteration at a time, so that greta can handle numerical errors as 'bad' proposals (see below).


an optional initials object (or list of initials objects of length chains) giving initial values for some or all of the variables in the model. These will be used as the starting point for sampling/optimisation.


the number of posterior samples to process at a time when tracing the parameters of interest; reduce this to reduce memory demands


a greta_mcmc_list object returned by mcmc or stashed_samples


named numeric values, giving initial values of some or all of the variables in the model (unnamed variables will be automatically initialised)


an optimiser object giving the optimisation algorithm and parameters See optimisers().


the maximum number of iterations before giving up


the numerical tolerance for the solution, the optimiser stops when the (absolute) difference in the joint density between successive iterations drops below this level


whether to account for Jacobian adjustments in the joint density. Set to FALSE (and do not use priors) for maximum likelihood estimates, or TRUE for maximum a posteriori estimates.


whether to return a list of analytically differentiated Hessian arrays for the parameters


For mcmc() if verbose = TRUE, the progress bar shows the number of iterations so far and the expected time to complete the phase of model fitting (warmup or sampling). Occasionally, a proposed set of parameters can cause numerical instability (I.e. the log density or its gradient is NA, Inf or -Inf); normally because the log joint density is so low that it can't be represented as a floating point number. When this happens, the progress bar will also display the proportion of proposals so far that were 'bad' (numerically unstable) and therefore rejected. Some numerical instability during the warmup phase is normal, but 'bad' samples during the sampling phase can lead to bias in your posterior sample. If you only have a few bad samples (<10\%), you can usually resolve this with a longer warmup period or by manually defining starting values to move the sampler into a more reasonable part of the parameter space. If you have more samples than that, it may be that your model is misspecified. You can often diagnose this by using calculate() to evaluate the values of greta arrays, given fixed values of model parameters, and checking the results are what you expect.

greta runs multiple chains simultaneously with a single sampler, vectorising all operations across the chains. E.g. a scalar addition in your model is computed as an elementwise vector addition (with vectors having length chains), a vector addition is computed as a matrix addition etc. TensorFlow is able to parallelise these operations, and this approach reduced computational overheads, so this is the most efficient of computing on multiple chains.

Multiple mcmc samplers (each of which can simultaneously run multiple chains) can also be run in parallel by setting the execution plan with the future package. Only plan(multisession) futures or plan(cluster) futures that don't use fork clusters are allowed, since forked processes conflict with TensorFlow's parallelism. Explicitly parallelising chains on a local machine with plan(multisession) will probably be slower than running multiple chains simultaneously in a single sampler (with plan(sequential), the default) because of the overhead required to start new sessions. However, plan(cluster) can be used to run chains on a cluster of machines on a local or remote network. See future::cluster() for details, and the future.batchtools package to set up plans on clusters with job schedulers.

If n_cores = NULL and mcmc samplers are being run sequentially, each sampler will be allowed to use all CPU cores (possibly to compute multiple chains sequentially). If samplers are being run in parallel with the future package, n_cores will be set so that ⁠n_cores * [future::nbrOfWorkers]⁠ is less than the number of CPU cores.

After carrying out mcmc on all the model parameters, mcmc() calculates the values of (i.e. traces) the parameters of interest for each of these samples, similarly to calculate(). Multiple posterior samples can be traced simultaneously, though this can require large amounts of memory for large models. As in calculate, the argument trace_batch_size can be modified to trade-off speed against memory usage.

If the sampler is aborted before finishing (and future parallelism isn't being used), the samples collected so far can be retrieved with stashed_samples(). Only samples from the sampling phase will be returned.

Samples returned by mcmc() and stashed_samples() can be added to with extra_samples(). This continues the chain from the last value of the previous chain and uses the same sampler and model as was used to generate the previous samples. It is not possible to change the sampler or extend the warmup period.

Because opt() acts on a list of greta arrays with possibly varying dimension, the par and hessian objects returned by opt() are named lists, rather than a vector (par) and a matrix (hessian), as returned by stats::optim(). Because greta arrays may not be vectors, the Hessians may not be matrices, but could be higher-dimensional arrays. To return a Hessian matrix covering multiple model parameters, you can construct your model so that all those parameters are in a vector, then split the vector up to define the model. The parameter vector can then be passed to model. See example.


mcmc, stashed_samples & extra_samples - a greta_mcmc_list object that can be analysed using functions from the coda package. This will contain mcmc samples of the greta arrays used to create model.

opt - a list containing the following named elements:

  • par a named list of the optimal values for the greta arrays specified in model

  • value the (unadjusted) negative log joint density of the model at the parameters 'par'

  • iterations the number of iterations taken by the optimiser

  • convergence an integer code, 0 indicates successful completion, 1 indicates the iteration limit max_iterations had been reached

  • hessian (if hessian = TRUE) a named list of hessian matrices/arrays for the parameters (w.r.t. value)


## Not run: 
# define a simple Bayesian model
x <- rnorm(10)
mu <- normal(0, 5)
sigma <- lognormal(1, 0.1)
distribution(x) <- normal(mu, sigma)
m <- model(mu, sigma)

# carry out mcmc on the model
draws <- mcmc(m, n_samples = 100)

# add some more samples
draws <- extra_samples(draws, 200)

#' # initial values can be passed for some or all model variables
draws <- mcmc(m, chains = 1, initial_values = initials(mu = -1))

# if there are multiple chains, a list of initial values should be passed,
# othewise the same initial values will be used for all chains
inits <- list(initials(sigma = 0.5), initials(sigma = 1))
draws <- mcmc(m, chains = 2, initial_values = inits)

# you can auto-generate a list of initials with something like this:
inits <- replicate(4,
  initials(mu = rnorm(1), sigma = runif(1)),
  simplify = FALSE
draws <- mcmc(m, chains = 4, initial_values = inits)

# or find the MAP estimate
opt_res <- opt(m)

# get the MLE of the normal variance
mu <- variable()
variance <- variable(lower = 0)
distribution(x) <- normal(mu, sqrt(variance))
m2 <- model(variance)

# adjust = FALSE skips the jacobian adjustments used in MAP estimation, to
# give the true maximum likelihood estimates
o <- opt(m2, adjust = FALSE)

# the MLE corresponds to the *unadjusted* sample variance, but differs
# from the sample variance
mean((x - mean(x))^2) # same
var(x) # different

# initial values can also be passed to optimisers:
o <- opt(m2, initial_values = initials(variance = 1))

# and you can return a list of the Hessians for each of these parameters
o <- opt(m2, hessian = TRUE)

# to get a hessian matrix across multiple greta arrays, you must first
# combine them and then split them up for use in the model (so that the
# combined vector is part of the model) and pass that vector to model:
params <- c(variable(), variable(lower = 0))
mu <- params[1]
variance <- params[2]
distribution(x) <- normal(mu, sqrt(variance))
m3 <- model(params)
o <- opt(m3, hessian = TRUE)

## End(Not run)

Install Python dependencies for greta


This is a helper function to install Python dependencies needed. This includes Tensorflow version 1.14.0, Tensorflow Probability 0.7.0, and numpy version 1.16.4. These Python modules will be installed into a virtual or conda environment, named "greta-env". Note that "virtualenv" is not available on Windows.


  method = c("auto", "virtualenv", "conda"),
  conda = "auto",
  timeout = 5,

  method = c("auto", "virtualenv", "conda"),
  conda = "auto",
  timeout = 5



Installation method ("virtualenv" or "conda")


The path to a conda executable. Use "auto" to allow reticulate to automatically find an appropriate conda binary. See Finding Conda for more details.


maximum time in minutes until the installation for each installation component times out and exits. Default is 5 minutes per installation component.


Optional arguments, reserved for future expansion.


This will automatically install Miniconda (a minimal version of the Anaconda scientific software management system), create a 'conda' environment for greta named 'greta-env' with required python and python package versions, and forcibly switch over to using that conda environment.

If you don't want to use conda or the "greta-env" conda environment, you can install these specific versions of tensorflow (version 1.14.0), and tensorflow-probability (version 0.7.0), and ensure that the python environment that is initialised in this R session has these versions installed. This is now always straightforward, so we recommend installing the python packages using install_greta_deps() for most users.


## Not run: 

## End(Not run)
## Not run: 
# to help troubleshoot your greta installation, this can help resolve some
# issues with installing greta dependencies

## End(Not run)

internal greta methods


A list of functions and R6 class objects that can be used to develop extensions to greta. Most users will not need to access these methods, and it is not recommended to use them directly in model code.


This help file lists the available internals, but they are not fully documented and are subject to change and deprecation without warning (though care will be taken not to break dependent packages on CRAN). For an overview of how greta works internally, see the technical details vignette. See for examples of R packages extending and building on greta.

Please get in contact via GitHub if you want to develop an extension to greta and need more details of how to use these internal functions.

You can use attach() to put a sublist in the search path. E.g. attach(.internals$nodes$constructors) will enable you to call op(), vble() and distrib() directly.


 .internals$greta_arrays$unknowns        # greta array print methods
 .internals$inference$progress_bar       # progress bar tools
                      samplers           # MCMC samplers
                      stash              # stashing MCMC samples
 .internals$nodes$constructors           # node creation wrappers
                  distribution_classes   # R6 distribution classes
                  mixture_classes        # R6 mixture distribution classes
                  node_classes           # R6 node classes
 .internals$tensors                      # functions on tensors
 .internals$utils$checks                 # checking function inputs
                  colours                # greta colour scheme
                  dummy_arrays           # mocking up extract/replace
                  misc                   # code simplification etc.
                  samplers               # mcmc helpers
 .internals$greta_stash                  # internal information storage

define joint distributions


joint combines univariate probability distributions together into a multivariate (and a priori independent between dimensions) joint distribution, either over a variable, or for fixed data.


joint(..., dim = NULL)



scalar variable greta arrays following probability distributions (see distributions()); the components of the joint distribution.


the dimensions of the greta array to be returned, either a scalar or a vector of positive integers. The final dimension of the greta array returned will be determined by the number of component distributions


The component probability distributions must all be either continuous or discrete, and must have the same dimensions.

This functionality is unlikely to be useful in most models, since the same result can usually be achieved by combining variables with separate distributions. It is included for situations where it is more convenient to consider these as a single distribution, e.g. for use with distribution or mixture.


## Not run: 
# an uncorrelated bivariate normal
x <- joint(normal(-3, 0.5), normal(3, 0.5))
m <- model(x)
plot(mcmc(m, n_samples = 500))

# joint distributions can be used to define densities over data
x <- cbind(rnorm(10, 2, 0.5), rbeta(10, 3, 3))
mu <- normal(0, 10)
sd <- normal(0, 3, truncation = c(0, Inf))
a <- normal(0, 3, truncation = c(0, Inf))
b <- normal(0, 3, truncation = c(0, Inf))
distribution(x) <- joint(normal(mu, sd), beta(a, b),
  dim = 10
m <- model(mu, sd, a, b)

## End(Not run)

mixtures of probability distributions


mixture combines other probability distributions into a single mixture distribution, either over a variable, or for fixed data.


mixture(..., weights, dim = NULL)



variable greta arrays following probability distributions (see distributions()); the component distributions in a mixture distribution.


a column vector or array of mixture weights, which must be positive, but need not sum to one. The first dimension must be the number of distributions, the remaining dimensions must either be 1 or match the distribution dimension.


the dimensions of the greta array to be returned, either a scalar or a vector of positive integers.


The weights are rescaled to sum to one along the first dimension, and are then used as the mixing weights of the distribution. I.e. the probability density is calculated as a weighted sum of the component probability distributions passed in via ⁠\dots⁠

The component probability distributions must all be either continuous or discrete, and must have the same dimensions.


## Not run: 
# a scalar variable following a strange bimodal distibution
weights <- uniform(0, 1, dim = 3)
a <- mixture(normal(-3, 0.5),
  normal(3, 0.5),
  normal(0, 3),
  weights = weights
m <- model(a)
plot(mcmc(m, n_samples = 500))

# simulate a mixture of poisson random variables and try to recover the
# parameters with a Bayesian model
x <- c(
  rpois(800, 3),
  rpois(200, 10)

weights <- uniform(0, 1, dim = 2)
rates <- normal(0, 10, truncation = c(0, Inf), dim = 2)
distribution(x) <- mixture(poisson(rates[1]),
  weights = weights
m <- model(rates)
draws_rates <- mcmc(m, n_samples = 500)

# check the mixing probabilities after fitting using calculate()
# (you could also do this within the model)
normalized_weights <- weights / sum(weights)
draws_weights <- calculate(normalized_weights, draws_rates)

# get the posterior means
summary(draws_rates)$statistics[, "Mean"]
summary(draws_weights)$statistics[, "Mean"]

# weights can also be an array, giving different mixing weights
# for each observation (first dimension must be number of components)
dim <- c(5, 4)
weights <- uniform(0, 1, dim = c(2, dim))
b <- mixture(normal(1, 1, dim = dim),
  normal(-1, 1, dim = dim),
  weights = weights

## End(Not run)

greta model objects


Create a greta_model object representing a statistical model (using model), and plot a graphical representation of the model. Statistical inference can be performed on greta_model objects with mcmc()


model(..., precision = c("double", "single"), compile = TRUE)

## S3 method for class 'greta_model'
print(x, ...)

## S3 method for class 'greta_model'
plot(x, y, colour = "#996bc7", ...)



for model: greta_array objects to be tracked by the model (i.e. those for which samples will be retained during mcmc). If not provided, all of the non-data greta_array objects defined in the calling environment will be tracked. For print and plot:further arguments passed to or from other methods (currently ignored).


the floating point precision to use when evaluating this model. Switching from "double" (the default) to "single" may decrease the computation time but increase the risk of numerical instability during sampling.


whether to apply XLA JIT compilation to the TensorFlow graph representing the model. This may slow down model definition, and speed up model evaluation.


a greta_model object


unused default argument


base colour used for plotting. Defaults to greta colours in violet.


model() takes greta arrays as arguments, and defines a statistical model by finding all of the other greta arrays on which they depend, or which depend on them. Further arguments to model can be used to configure the TensorFlow graph representing the model, to tweak performance.

The plot method produces a visual representation of the defined model. It uses the DiagrammeR package, which must be installed first. Here's a key to the plots: plotlegend.png


model - a greta_model object.

plot - a DiagrammeR::grViz() object, with the DiagrammeR::dgr_graph() object used to create it as an attribute "dgr_graph".


## Not run: 

# define a simple model
mu <- variable()
sigma <- normal(0, 3, truncation = c(0, Inf))
x <- rnorm(10)
distribution(x) <- normal(mu, sigma)

m <- model(mu, sigma)


## End(Not run)

arithmetic, logical and relational operators for greta arrays


This is a list of currently implemented arithmetic, logical and relational operators to combine greta arrays into probabilistic models. Also see functions and transforms.


greta's operators are used just like R's the standard arithmetic, logical and relational operators, but they return other greta arrays. Since the operations are only carried during sampling, the greta array objects have unknown values.


 # arithmetic operators
 x + y
 x - y
 x * y
 x / y
 x ^ y
 x %% y
 x %/% y
 x %*% y

 # logical operators
 x & y
 x | y

 # relational operators
 x < y
 x > y
 x <= y
 x >= y
 x == y
 x != y


## Not run: 

x <- as_data(-1:12)

# arithmetic
a <- x + 1
b <- 2 * x + 3
c <- x %% 2
d <- x %/% 5

# logical
e <- (x > 1) | (x < 1)
f <- e & (x < 2)
g <- !f

# relational
h <- x < 1
i <- (-x) >= x
j <- h == x

## End(Not run)

optimisation methods


Functions to set up optimisers (which find parameters that maximise the joint density of a model) and change their tuning parameters, for use in opt(). For details of the algorithms and how to tune them, see the SciPy optimiser docs or the TensorFlow optimiser docs.







l_bfgs_b(maxcor = 10, maxls = 20)

tnc(max_cg_it = -1, stepmx = 0, rescale = -1)

cobyla(rhobeg = 1)


gradient_descent(learning_rate = 0.01)

adadelta(learning_rate = 0.001, rho = 1, epsilon = 1e-08)

adagrad(learning_rate = 0.8, initial_accumulator_value = 0.1)

  learning_rate = 0.8,
  global_step = 1L,
  initial_gradient_squared_accumulator_value = 0.1,
  l1_regularization_strength = 0,
  l2_regularization_strength = 0

momentum(learning_rate = 0.001, momentum = 0.9, use_nesterov = TRUE)

adam(learning_rate = 0.1, beta1 = 0.9, beta2 = 0.999, epsilon = 1e-08)

  learning_rate = 1,
  learning_rate_power = -0.5,
  initial_accumulator_value = 0.1,
  l1_regularization_strength = 0,
  l2_regularization_strength = 0

  learning_rate = 0.01,
  l1_regularization_strength = 0,
  l2_regularization_strength = 0

  learning_rate = 1,
  initial_accumulator_value = 0.1,
  l1_regularization_strength = 0,
  l2_regularization_strength = 0

rms_prop(learning_rate = 0.1, decay = 0.9, momentum = 0, epsilon = 1e-10)



maximum number of 'variable metric corrections' used to define the approximation to the hessian matrix


maximum number of line search steps per iteration


maximum number of hessian * vector evaluations per iteration


maximum step for the line search


log10 scaling factor used to trigger rescaling of objective


reasonable initial changes to the variables


the size of steps (in parameter space) towards the optimal value


the decay rate


a small constant used to condition gradient updates


initial value of the 'accumulator' used to tune the algorithm


the current training step number


initial value of the accumulators used to tune the algorithm


L1 regularisation coefficient (must be 0 or greater)


L2 regularisation coefficient (must be 0 or greater)


the momentum of the algorithm


whether to use Nesterov momentum


exponential decay rate for the 1st moment estimates


exponential decay rate for the 2nd moment estimates


power on the learning rate, must be 0 or less


discounting factor for the gradient


The optimisers powell(), cg(), newton_cg(), l_bfgs_b(), tnc(), cobyla(), and slsqp() are deprecated. They will be removed in greta 0.4.0, since they will no longer be available in TensorFlow 2.0, on which that version of greta will depend.

The cobyla() does not provide information about the number of iterations nor convergence, so these elements of the output are set to NA


an optimiser object that can be passed to opt().


## Not run: 
# use optimisation to find the mean and sd of some data
x <- rnorm(100, -2, 1.2)
mu <- variable()
sd <- variable(lower = 0)
distribution(x) <- normal(mu, sd)
m <- model(mu, sd)

# configure optimisers & parameters via 'optimiser' argument to opt
opt_res <- opt(m, optimiser = bfgs())

# compare results with the analytic solution
c(mean(x), sd(x))

## End(Not run)

Functions overloaded by greta


greta provides a wide range of methods to apply common R functions and operations to greta_array objects. A few of these functions and operators are not associated with a class system, so they are overloaded here. This should not affect normal use of these functions, but they need to be documented to satisfy CRAN's check.


x %*% y

chol2inv(x, size = NCOL(x), LINPACK = FALSE)



colMeans(x, na.rm = FALSE, dims = 1L)

rowMeans(x, na.rm = FALSE, dims = 1L)

colSums(x, na.rm = FALSE, dims = 1L)

rowSums(x, na.rm = FALSE, dims = 1L)

sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...)

backsolve(r, x, k = ncol(r), upper.tri = TRUE, transpose = FALSE)

forwardsolve(l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)

apply(X, MARGIN, FUN, ...)

tapply(X, INDEX, FUN, ...)

eigen(x, symmetric, only.values, EISPACK)

rdist(x1, x2 = NULL, compact = FALSE)

  along = N,
  rev.along = NULL,
  new.names = NULL,
  force.array = TRUE,
  make.names = use.anon.names,
  use.anon.names = FALSE,
  use.first.dimnames = FALSE,
  hier.names = FALSE,
  use.dnns = FALSE

diag(x = 1, nrow, ncol)


x, y, size, LINPACK, V, na.rm, dims, MARGIN, STATS, FUN, check.margin, ..., r, k, upper.tri, transpose, l, X, INDEX, symmetric, only.values, EISPACK, x1, x2, compact, along, rev.along, new.names, force.array, make.names, use.anon.names, use.first.dimnames, hier.names, use.dnns, nrow, ncol

arguments as in original documentation


Note that, since R 3.1, the LINPACK argument is defunct and silently ignored. The argument is only included for compatibility with the base functions that call it.

To find the original help file for these overloaded functions, search for the function, e.g., ?cov2cor and select the non-greta function.

Helpers to remove, and reinstall python environments and miniconda


This can be useful when debugging greta installation to get to "clean slate". There are four functions:



reinstall_greta_env(timeout = 5)


reinstall_miniconda(timeout = 5)



time in minutes to wait until timeout (default is 5 minutes)


  • remove_greta_env() removes the 'greta-env' conda environment

  • remove_miniconda() removes miniconda installation

  • reinstall_greta_env() remove 'greta-env' and reinstall it using greta_create_conda_env() (which is used internally).

  • reinstall_miniconda() removes miniconda and reinstalls it using greta_install_miniconda() (which is used internally)




## Not run: 

## End(Not run)

MCMC samplers


Functions to set up MCMC samplers and change the starting values of their parameters, for use in mcmc().


hmc(Lmin = 5, Lmax = 10, epsilon = 0.1, diag_sd = 1)

rwmh(proposal = c("normal", "uniform"), epsilon = 0.1, diag_sd = 1)

slice(max_doublings = 5)



minimum number of leapfrog steps (positive integer, Lmin > Lmax)


maximum number of leapfrog steps (positive integer, Lmax > Lmin)


leapfrog stepsize hyperparameter (positive, will be tuned)


estimate of the posterior marginal standard deviations (positive, will be tuned).


the probability distribution used to generate proposal states


the maximum number of iterations of the 'doubling' algorithm used to adapt the size of the slice


During the warmup iterations of mcmc, some of these sampler parameters will be tuned to improve the efficiency of the sampler, so the values provided here are used as starting values.

For hmc(), the number of leapfrog steps at each iteration is selected uniformly at random from between Lmin and Lmax. diag_sd is used to rescale the parameter space to make it more uniform, and make sampling more efficient.

rwmh() creates a random walk Metropolis-Hastings sampler; a a gradient-free sampling algorithm. The algorithm involves a proposal generating step proposal_state = current_state + perturb by a random perturbation, followed by Metropolis-Hastings accept/reject step. The class is implemented for uniform and normal proposals.

slice() implements a multivariate slice sampling algorithm. Currently this algorithm can only be used with single-precision models (set using the precision argument to model()). The parameter max_doublings is not tuned during warmup.


a sampler object that can be passed to mcmc().

Simulate Responses From greta_model Object


Simulate values of all named greta arrays associated with a greta model from the model priors, including the response variable.


## S3 method for class 'greta_model'
simulate(object, nsim = 1, seed = NULL, precision = c("double", "single"), ...)



a greta_model() object


positive integer scalar - the number of responses to simulate


an optional seed to be used in set.seed immediately before the simulation so as to generate a reproducible sample


the floating point precision to use when calculating values.


optional additional arguments, none are used at present


This is essentially a wrapper around calculate() that finds all relevant greta arrays. See that function for more functionality, including simulation conditional on fixed values or posterior samples.

To simulate values of the response variable, it must be both a named object (in the calling environment) and be a greta array. If you don't see it showing up in the output, you may need to use as_data to convert it to a greta array before defining the model.


A named list of vectors, matrices or arrays containing independent samples of the greta arrays associated with the model. The number of samples will be prepended as the first dimension of the greta array, so that a vector of samples is returned for each scalar greta array, and a matrix is returned for each vector greta array, etc.


## Not run: 
# build a greta model
n <- 10
y <- rnorm(n)
y <- as_data(y)

sd <- lognormal(1, 2)
mu <- normal(0, 1, dim = n)
distribution(y) <- normal(mu, sd)
m <- model(mu, sd)

# simulate one random draw of y, mu and sd from the model prior:
sims <- simulate(m)

# 100 simulations of y, mu and sd
sims <- simulate(m, nsim = 100)

## End(Not run)
# nolint start

create data greta arrays


These structures can be used to set up more complex models. For example, scalar parameters can be embedded in a greta array by first creating a greta array with zeros() or ones(), and then embedding the parameter value using greta's replacement syntax.




greta_array(data = 0, dim = length(data))



dimensions of the greta arrays to create


a vector giving data to fill the greta array. Other object types are coerced by as.vector().


an integer vector giving the dimensions for the greta array to be created.


greta_array is a convenience function to create an R array with array() and then coerce it to a greta array. I.e. when passed something that can be coerced to a numeric array, it is equivalent to as_data(array(data, dim)).

If data is a greta array and dim is different than dim(data), a reshaped greta array is returned. This is equivalent to: dim(data) <- dim.


a greta array object


## Not run: 

# a 3 row, 4 column greta array of 0s
z <- zeros(3, 4)

# a 3x3x3 greta array of 1s
z <- ones(3, 3, 3)

# a 2x4 greta array filled with pi
z <- greta_array(pi, dim = c(2, 4))

# a 3x3x3 greta array filled with 1, 2, and 3
z <- greta_array(1:3, dim = c(3, 3, 3))

## End(Not run)

transformation functions for greta arrays


transformations for greta arrays, which may also be used as inverse link functions. Also see operators and functions.










a real-valued (i.e. values ranging from -Inf to Inf) greta array to transform to a constrained value


greta does not allow you to state the transformation/link on the left hand side of an assignment, as is common in the BUGS and STAN modelling languages. That's because the same syntax has a very different meaning in R, and can only be applied to objects that are already in existence. The inverse forms of the common link functions (prefixed with an 'i') can be used instead.

The log1pe inverse link function is equivalent to log(1 + exp(x)), yielding a positive transformed parameter. Unlike the log transformation, this transformation is approximately linear for x > 1. i.e. when x>1x > 1, yy is approximately xx

imultilogit expects an n-by-m greta array, and returns an n-by-(m+1) greta array of positive reals whose rows sum to one. This is equivalent adding a final column of 0s and then running the softmax function widely used in machine learning.


## Not run: 

x1 <- normal(1, 3, dim = 10)

# transformation to the unit interval
p1 <- iprobit(x1)
p2 <- ilogit(x1)
p3 <- icloglog(x1)
p4 <- icauchit(x1)

# and to positive reals
y <- log1pe(x1)

# transform from 10x3 to 10x4, where rows are a complete set of
# probabilities
x2 <- normal(1, 3, dim = c(10, 3))
z <- imultilogit(x2)

## End(Not run)

create greta variables


variable() creates greta arrays representing unknown parameters, to be learned during model fitting. These parameters are not associated with a probability distribution. To create a variable greta array following a specific probability distribution, see distributions().


variable(lower = -Inf, upper = Inf, dim = NULL)

cholesky_variable(dim, correlation = FALSE)




lower, upper

optional limits to variables. These must be specified as numerics, they cannot be greta arrays (though see details for a workaround). They can be set to -Inf (lower) or Inf (upper), though lower must always be less than upper.


the dimensions of the greta array to be returned, either a scalar or a vector of positive integers. See details.


whether to return a cholesky factor corresponding to a correlation matrix (diagonal elements equalling 1, off-diagonal elements between -1 and 1).


lower and upper must be fixed, they cannot be greta arrays. This ensures these values can always be transformed to a continuous scale to run the samplers efficiently. However, a variable parameter with dynamic limits can always be created by first defining a variable constrained between 0 and 1, and then transforming it to the required scale. See below for an example.

The constraints in simplex_variable() and ordered_variable() operate on the final dimension, which must have more than 1 element. Passing in a scalar value for dim therefore results in a row-vector.


## Not run: 

# a scalar variable
a <- variable()

# a positive length-three variable
b <- variable(lower = 0, dim = 3)

# a 2x2x2 variable bounded between 0 and 1
c <- variable(lower = 0, upper = 1, dim = c(2, 2, 2))

# create a variable, with lower and upper defined by greta arrays
min <- as_data(iris$Sepal.Length)
max <- min^2
d <- min + variable(0, 1, dim = nrow(iris)) * (max - min)

## End(Not run)
# 4x4 cholesky factor variables for covariance and correlation matrices
e_cov <- cholesky_variable(dim = 4)
e_correl <- cholesky_variable(dim = 4, correlation = TRUE)

# these can be converted to symmetic matrices with chol2symm
# (equivalent to t(e_cov) %*% e_cov, but more efficient)
cov <- chol2symm(e_cov)
correl <- chol2symm(e_correl)
# a 4D simplex (sums to 1, all values positive)
f <- simplex_variable(4)

# a 4D simplex on the final dimension
g <- simplex_variable(dim = c(2, 3, 4))
# a 2D variable with each element higher than the one in the cell to the left
h <- ordered_variable(dim = c(3, 4))

# more constraints can be added with monotonic transformations, e.g. an
# ordered positive variable
i <- exp(ordered_variable(5))