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:  20240709 05:05:32 UTC 
Source:  https://github.com/gretadev/greta 
define an object in an R session as a data greta array for use as data in a greta model.
as_data(x)
as_data(x)
x 
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)
## 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 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.
calculate( ..., values = list(), nsim = NULL, seed = NULL, precision = c("double", "single"), trace_batch_size = 100 )
calculate( ..., 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 
values 
a named list giving temporary values of the greta arrays with
which 
nsim 
an optional positive integer scalar for the number of responses to simulate if stochastic greta arrays are present in the model  see Details. 
seed 
an optional seed to be used in set.seed immediately before the simulation so as to generate a reproducible sample 
precision 
the floating point precision to use when calculating values. 
trace_batch_size 
the number of posterior samples to process at a time
when 
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) calculate(y) # 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), max(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)
## 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) calculate(y) # 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), max(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)
Evaluate t(x) \%*\% x
efficiently, where x
is the
(uppertriangular) Cholesky factor of a symmetric, positive definite square
matrix. I.e. it is the inverse of chol
.
chol2symm(x)
chol2symm(x)
x 
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)
# 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)
distribution
defines probability distributions over
observed data, e.g. to set a model likelihood.
distribution(greta_array) < value distribution(greta_array)
distribution(greta_array) < value distribution(greta_array)
greta_array 
a data greta array. For the assignment method it must not already have a probability distribution assigned 
value 
a greta array with a distribution (see

The extract method returns the greta array if it has a distribution,
or NULL
if it doesn't. It has no real usecase, 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 distribution(y) ## End(Not run)
## 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 distribution(y) ## End(Not run)
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)
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 
dim 
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, 
truncation 
a lengthtwo vector giving values between which to truncate
the distribution, similarly to the 
prob 
probability parameter ( 
size , m , n , k

positive integer parameter 
Sigma 
positive definite variancecovariance matrix parameter 
n_realisations 
the number of independent realisation of a multivariate distribution 
dimension 
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 lengththree variable, with each element following a standard normal # distribution alpha < normal(0, 1, dim = 3) # a lengththree 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) dim(theta) # a Wishart variable with the same covariance parameter theta < wishart(df = 5, Sigma = Sig) ## End(Not run)
## Not run: # a uniform parameter constrained to be between 0 and 1 phi < uniform(min = 0, max = 1) # a lengththree variable, with each element following a standard normal # distribution alpha < normal(0, 1, dim = 3) # a lengththree 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) dim(theta) # a Wishart variable with the same covariance parameter theta < wishart(df = 5, Sigma = Sig) ## End(Not run)
Generic methods to extract and replace elements of greta arrays, or to combine greta arrays.
x 
a greta array 
i , j

indices specifying elements to extract or replace 
n 
a single integer, as in 
nrow , ncol

optional dimensions for the resulting greta array when x is not a matrix. 
value 
for 
... 
either further indices specifying elements to extract or replace
( 
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 twodimensional greta array, but it cannot be used to create a
matrixlike greta array from a scalar or vectorlike 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] 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 cbind(...) rbind(...) abind(...) c(..., recursive = FALSE) rep(x, times, ..., recursive = FALSE) # get and set dimensions length(x) dim(x) 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)
## 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)
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 twodimensional 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 log(x) exp(x) log1p(x) expm1(x) # miscellaneous mathematics abs(x) mean(x) sqrt(x) sign(x) # rounding of numbers ceiling(x) floor(x) round(x, digits = 0) # trigonometry cos(x) sin(x) tan(x) acos(x) asin(x) atan(x) cosh(x) sinh(x) tanh(x) acosh(x) asinh(x) atanh(x) cospi(x) sinpi(x) tanpi(x) # special mathematical functions lgamma(x) digamma(x) trigamma(x) choose(n, k) lchoose(n, k) # matrix operations t(x) chol(x, ...) chol2inv(x, ...) cov2cor(V) 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 cumsum(x) cumprod(x) cummax(x) cummin(x) # 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)
## 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 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)
Authors:
Nick Golding [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]
Useful links:
Report bugs at https://github.com/gretadev/greta/issues
## 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)
## 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)
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.
greta_notes_install_miniconda_output() greta_notes_install_miniconda_error() greta_notes_conda_create_output() greta_notes_conda_create_error() greta_notes_conda_install_output() greta_notes_conda_install_error() greta_notes_tf_num_error()
greta_notes_install_miniconda_output() greta_notes_install_miniconda_error() greta_notes_conda_create_output() greta_notes_conda_create_error() greta_notes_conda_install_output() greta_notes_conda_install_error() greta_notes_tf_num_error()
## Not run: greta_notes_install_miniconda() greta_notes_conda_create() greta_notes_conda_install() greta_notes_tf_num_error() greta_notes_tf_error() ## End(Not run)
## Not run: greta_notes_install_miniconda() greta_notes_conda_create() greta_notes_conda_install() greta_notes_tf_num_error() greta_notes_tf_error() ## End(Not run)
This checks if Python, Tensorflow, Tensorflow Probability, and the greta conda environment are available, and also loads and initialises python
greta_sitrep()
greta_sitrep()
Message if greta is ready to use
## Not run: greta_sitrep() ## End(Not run)
## Not run: greta_sitrep() ## End(Not run)
Carry out statistical inference on greta models by MCMC or likelihood/posterior optimisation.
mcmc( model, 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 ) stashed_samples() extra_samples( draws, n_samples = 1000, thin = 1, n_cores = NULL, verbose = TRUE, pb_update = 50, one_by_one = FALSE, trace_batch_size = 100 ) initials(...) opt( model, optimiser = bfgs(), max_iterations = 100, tolerance = 1e06, initial_values = initials(), adjust = TRUE, hessian = FALSE )
mcmc( model, 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 ) stashed_samples() extra_samples( draws, n_samples = 1000, thin = 1, n_cores = NULL, verbose = TRUE, pb_update = 50, one_by_one = FALSE, trace_batch_size = 100 ) initials(...) opt( model, optimiser = bfgs(), max_iterations = 100, tolerance = 1e06, initial_values = initials(), adjust = TRUE, hessian = FALSE )
model 
greta_model object 
sampler 
sampler used to draw values in MCMC. See

n_samples 
number of MCMC samples to draw per chain (after any warmup, but before thinning) 
thin 
MCMC thinning rate; every 
warmup 
number of samples to spend warming up the mcmc sampler (moving chains toward the highest density area and tuning sampler hyperparameters). 
chains 
number of MCMC chains to run 
n_cores 
the maximum number of CPU cores used by each sampler (see details). 
verbose 
whether to print progress information to the console 
pb_update 
how regularly to update the progress bar (in iterations).
If 
one_by_one 
whether to run TensorFlow MCMC code one iteration at a time, so that greta can handle numerical errors as 'bad' proposals (see below). 
initial_values 
an optional 
trace_batch_size 
the number of posterior samples to process at a time when tracing the parameters of interest; reduce this to reduce memory demands 
draws 
a greta_mcmc_list object returned by 
... 
named numeric values, giving initial values of some or all of the variables in the model (unnamed variables will be automatically initialised) 
optimiser 
an 
max_iterations 
the maximum number of iterations before giving up 
tolerance 
the numerical tolerance for the solution, the optimiser stops when the (absolute) difference in the joint density between successive iterations drops below this level 
adjust 
whether to account for Jacobian adjustments in the joint
density. Set to 
hessian 
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 tradeoff 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 higherdimensional 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 autogenerate 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 o$par 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) o$hessian # 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) o$hessian ## End(Not run)
## 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 autogenerate 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 o$par 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) o$hessian # 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) o$hessian ## End(Not run)
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 "gretaenv". Note that "virtualenv" is not available on Windows.
install_greta_deps( method = c("auto", "virtualenv", "conda"), conda = "auto", timeout = 5, ... ) reinstall_greta_deps( method = c("auto", "virtualenv", "conda"), conda = "auto", timeout = 5 )
install_greta_deps( method = c("auto", "virtualenv", "conda"), conda = "auto", timeout = 5, ... ) reinstall_greta_deps( method = c("auto", "virtualenv", "conda"), conda = "auto", timeout = 5 )
method 
Installation method ("virtualenv" or "conda") 
conda 
The path to a 
timeout 
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 'gretaenv' 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 "gretaenv" conda environment, you
can install these specific versions of tensorflow (version 1.14.0), and
tensorflowprobability (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: install_greta_deps() ## End(Not run) ## Not run: # to help troubleshoot your greta installation, this can help resolve some # issues with installing greta dependencies reinstall_greta_deps() ## End(Not run)
## Not run: install_greta_deps() ## End(Not run) ## Not run: # to help troubleshoot your greta installation, this can help resolve some # issues with installing greta dependencies reinstall_greta_deps() ## End(Not run)
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 https://github.com/gretadev 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
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)
joint(..., dim = NULL)
... 
scalar variable greta arrays following probability distributions
(see 
dim 
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) plot(mcmc(m)) ## End(Not run)
## 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) plot(mcmc(m)) ## End(Not run)
mixture
combines other probability distributions into a
single mixture distribution, either over a variable, or for fixed data.
mixture(..., weights, dim = NULL)
mixture(..., weights, dim = NULL)
... 
variable greta arrays following probability distributions (see

weights 
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. 
dim 
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]), poisson(rates[2]), 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)
## 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]), poisson(rates[2]), 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)
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", ...)
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 
precision 
the floating point precision to use when evaluating this
model. Switching from 
compile 
whether to apply XLA JIT compilation to the TensorFlow graph representing the model. This may slow down model definition, and speed up model evaluation. 
x 
a 
y 
unused default argument 
colour 
base colour used for plotting. Defaults to 
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:
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) plot(m) ## End(Not run)
## 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) plot(m) ## End(Not run)
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 x + y x  y x * y x / y x ^ y x %% y x %/% y x %*% y # logical operators !x 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)
## 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)
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.
nelder_mead() powell() cg() bfgs() newton_cg() l_bfgs_b(maxcor = 10, maxls = 20) tnc(max_cg_it = 1, stepmx = 0, rescale = 1) cobyla(rhobeg = 1) slsqp() gradient_descent(learning_rate = 0.01) adadelta(learning_rate = 0.001, rho = 1, epsilon = 1e08) adagrad(learning_rate = 0.8, initial_accumulator_value = 0.1) adagrad_da( 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 = 1e08) ftrl( learning_rate = 1, learning_rate_power = 0.5, initial_accumulator_value = 0.1, l1_regularization_strength = 0, l2_regularization_strength = 0 ) proximal_gradient_descent( learning_rate = 0.01, l1_regularization_strength = 0, l2_regularization_strength = 0 ) proximal_adagrad( 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 = 1e10)
nelder_mead() powell() cg() bfgs() newton_cg() l_bfgs_b(maxcor = 10, maxls = 20) tnc(max_cg_it = 1, stepmx = 0, rescale = 1) cobyla(rhobeg = 1) slsqp() gradient_descent(learning_rate = 0.01) adadelta(learning_rate = 0.001, rho = 1, epsilon = 1e08) adagrad(learning_rate = 0.8, initial_accumulator_value = 0.1) adagrad_da( 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 = 1e08) ftrl( learning_rate = 1, learning_rate_power = 0.5, initial_accumulator_value = 0.1, l1_regularization_strength = 0, l2_regularization_strength = 0 ) proximal_gradient_descent( learning_rate = 0.01, l1_regularization_strength = 0, l2_regularization_strength = 0 ) proximal_adagrad( 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 = 1e10)
maxcor 
maximum number of 'variable metric corrections' used to define the approximation to the hessian matrix 
maxls 
maximum number of line search steps per iteration 
max_cg_it 
maximum number of hessian * vector evaluations per iteration 
stepmx 
maximum step for the line search 
rescale 
log10 scaling factor used to trigger rescaling of objective 
rhobeg 
reasonable initial changes to the variables 
learning_rate 
the size of steps (in parameter space) towards the optimal value 
rho 
the decay rate 
epsilon 
a small constant used to condition gradient updates 
initial_accumulator_value 
initial value of the 'accumulator' used to tune the algorithm 
global_step 
the current training step number 
initial_gradient_squared_accumulator_value 
initial value of the accumulators used to tune the algorithm 
l1_regularization_strength 
L1 regularisation coefficient (must be 0 or greater) 
l2_regularization_strength 
L2 regularisation coefficient (must be 0 or greater) 
momentum 
the momentum of the algorithm 
use_nesterov 
whether to use Nesterov momentum 
beta1 
exponential decay rate for the 1st moment estimates 
beta2 
exponential decay rate for the 2nd moment estimates 
learning_rate_power 
power on the learning rate, must be 0 or less 
decay 
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 opt_res$par c(mean(x), sd(x)) ## End(Not run)
## 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 opt_res$par c(mean(x), sd(x)) ## End(Not run)
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) cov2cor(V) identity(x) 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) abind( ..., 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 chol2inv(x, size = NCOL(x), LINPACK = FALSE) cov2cor(V) identity(x) 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) abind( ..., 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 nongreta function.
This can be useful when debugging greta installation to get to "clean slate". There are four functions:
remove_greta_env() reinstall_greta_env(timeout = 5) remove_miniconda() reinstall_miniconda(timeout = 5)
remove_greta_env() reinstall_greta_env(timeout = 5) remove_miniconda() reinstall_miniconda(timeout = 5)
timeout 
time in minutes to wait until timeout (default is 5 minutes) 
remove_greta_env()
removes the 'gretaenv' conda environment
remove_miniconda()
removes miniconda installation
reinstall_greta_env()
remove 'gretaenv' 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)
invisible
## Not run: remove_greta_env() remove_miniconda() reinstall_greta_env() reinstall_miniconda() ## End(Not run)
## Not run: remove_greta_env() remove_miniconda() reinstall_greta_env() reinstall_miniconda() ## End(Not run)
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)
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)
Lmin 
minimum number of leapfrog steps (positive integer, Lmin > Lmax) 
Lmax 
maximum number of leapfrog steps (positive integer, Lmax > Lmin) 
epsilon 
leapfrog stepsize hyperparameter (positive, will be tuned) 
diag_sd 
estimate of the posterior marginal standard deviations (positive, will be tuned). 
proposal 
the probability distribution used to generate proposal states 
max_doublings 
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 MetropolisHastings sampler; a
a gradientfree sampling algorithm. The algorithm involves a proposal
generating step proposal_state = current_state + perturb
by a random
perturbation, followed by MetropolisHastings 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 singleprecision 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()
.
greta_model
ObjectSimulate 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"), ...)
## S3 method for class 'greta_model' simulate(object, nsim = 1, seed = NULL, precision = c("double", "single"), ...)
object 
a 
nsim 
positive integer scalar  the number of responses to simulate 
seed 
an optional seed to be used in set.seed immediately before the simulation so as to generate a reproducible sample 
precision 
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) library(greta) 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
## Not run: # 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) # 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
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.
zeros(...) ones(...) greta_array(data = 0, dim = length(data))
zeros(...) ones(...) greta_array(data = 0, dim = length(data))
... 
dimensions of the greta arrays to create 
data 
a vector giving data to fill the greta array. Other object types
are coerced by 
dim 
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)
## 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)
transformations for greta arrays, which may also be used as inverse link functions. Also see operators and functions.
iprobit(x) ilogit(x) icloglog(x) icauchit(x) log1pe(x) imultilogit(x)
iprobit(x) ilogit(x) icloglog(x) icauchit(x) log1pe(x) imultilogit(x)
x 
a realvalued (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 > 1$
, $y$
is approximately $x$
imultilogit
expects an nbym greta array, and returns an nby(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)
## 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)
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) simplex_variable(dim) ordered_variable(dim)
variable(lower = Inf, upper = Inf, dim = NULL) cholesky_variable(dim, correlation = FALSE) simplex_variable(dim) ordered_variable(dim)
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 
dim 
the dimensions of the greta array to be returned, either a scalar or a vector of positive integers. See details. 
correlation 
whether to return a cholesky factor corresponding to a correlation matrix (diagonal elements equalling 1, offdiagonal 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 rowvector.
## Not run: # a scalar variable a < variable() # a positive lengththree 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))
## Not run: # a scalar variable a < variable() # a positive lengththree 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))