### Variables and priors The second section of the script creates three greta arrays to represent the parameters in our model: ```{r variables} int <- normal(0, 1) coef <- normal(0, 3) sd <- student(3, 0, 1, truncation = c(0, Inf)) ``` Each of these is a variable greta array, and each is assumed *a priori* (before fitting to data) to follow a different probability distribution. In other words, these are prior distributions over variables, which we need to specify to make this a full Bayesian analysis. Before going through how to specify variables with probability distributions, it will be clearer to first demonstrate the alternative; variables without probability distributions. #### variables without probability distributions If we were carrying out a frequentist analysis of this model, we could create variable greta arrays (values we want to learn) without probability distributions using the `variable()` function. E.g. in a frequentist version of the model we could create `int` with: ```{r int_variable} (int <- variable()) ``` `variable()` has three arguments. The first two arguments determine the constraints on this parameter; we left them at their default setting of `lower = -Inf, upper = Inf` meaning the variables can take any value on the real number line. The third argument gives the dimensions of the greta array to return, in this case we left it at its default value of 1x1 (a scalar). We can create a variable constrained between two values by specifying `lower` and `upper`. So we could have created the positive variable `sd` (the standard deviation of the likelihood) with: ```{r positive_variable} (sd <- variable(lower = 0)) ``` If we had instead wanted a 2x3 array of positive variables we could have created it like this: ```{r matrix_variable} variable(lower = 0, dim = c(2, 3)) ``` #### variables with probability distributions In our example script, when we created the variables `int`, `coef` and `sd`, we simultaneously stated the prior distributions for them using some of greta's probability distribution functions. You can see a list of the currently available distributions in the [`?greta::distributions`](https://greta-stats.org/reference/distributions.html) helpfile. Each of these distribution functions takes as arguments the distribution's parameters (either as numbers or other greta arrays), as well as the dimensions of the resulting greta array. As before, we left the dimensions arguments at their default value to create scalar greta arrays. Both `int` and `coef` were given zero-mean normal distributions, which are a common choice of prior for unconstrained variables in Bayesian analyses. For the strictly positive parameter `sd`, we chose a slightly unconventional option, a positive-truncated (non-standard) student's t distribution, which we create using greta's built-in support for truncated distributions. #### variables with truncated distributions Some of greta's probability distributions (those that are continuous and univariate) can be specified as truncated distributions. By modifying the `truncation` argument, we can state that the resulting distribution should be truncated between the two truncation bounds. So to create a standard normal distribution truncated between -1 and 1 we can do: ```{r truncated1} (z <- normal(0, 1, truncation = c(-1, 1))) ``` greta will account for this truncation when calculating the density of this distribution; rescaling it to be a valid probability distribution. We can only truncate to within the support of the distribution; e.g. greta will throw an error if we try to truncate a lognormal distribution (which must be positive) to have a lower bound of -1.
### Operations The third section of the script uses mathematical operators to combine some of our parameters and data, to calculate the predicted sepal lengths, for a given parameter set: ```{r linear_predictor} mean <- int + coef * x ``` Because `int` and `coef` are both scalars, the resulting greta array `mean` has the same dimensions as `x`; a column vector with 150 elements: ```{r mean} dim(mean) head(mean) ``` greta arrays can be manipulated using R's standard arithmetic, logical and relational operators; including `+`, `*` and many others. The [`?greta::operators`](https://greta-stats.org/reference/operators.html) help file lists the operators that are implemented for greta arrays. You can also use a lot of common R functions for numeric data, such as `sum()`, `log()` and others. the available functions are listed in the [`?greta::functions`](https://greta-stats.org/reference/functions.html) helpfile. All of these mathematical manipulations of greta arrays produce 'operation' greta arrays. #### Extract and replace We can use R's extract and replace syntax (using `[`) on greta arrays, just like with R's vectors, matrices and arrays. E.g. to extract elements from `mean` we can do: ```{r extract} mean[1:3] ``` We can assign values from one greta array to another too. For example, if we wanted to create a matrix that has random normal variables in the first column, but zeros everywhere else, we could do: ```{r replace} z <- zeros(4, 3) z[, 1] <- normal(0, 1, dim = 4) z ``` R's subset operator `[` has an argument `drop`, which determines whether to reduce the number of dimensions of a array or matrix when the object has zero elements in that dimension. By default, `drop = TRUE` for R objects, so matrices are automatically converted into vectors (which have dimension `NULL`) if you take out a single column: ```{r drop} z <- matrix(1, nrow = 2, ncol = 2) dim(z[, 1]) dim(z[, 1, drop = FALSE]) ``` greta arrays must always have two dimensions, so greta always acts as though `drop = FALSE`: ```{r drop_greta} z_greta <- as_data(z) dim(z_greta[, 1]) ``` #### Functions We can write our own functions for greta arrays, using the existing operators and functions. For example, we could define the [inverse hyperbolic tangent function](https://en.wikipedia.org/wiki/Inverse_hyperbolic_function#Inverse_hyperbolic_tangent) for greta arrays like this: ```{r function1} atanh <- function (z) (log(1 + z) - log(1 - z)) / 2 atanh(z_greta) ```
### Likelihood So far, we have created greta arrays representing the variables in our model (with prior distributions) and created other greta arrays from them and some fixed, *independent* data. To perform statistical inference on this model, we also need to link it to some observed *dependent* data. By comparing the sepal lengths predicted under different parameter values with the observed sepal lengths, we can estimate the most plausible values of those parameters. We do that by defining a likelihood for the observed sepal length data `y`. By defining a likelihood over observed data, we are stating that these observed data are actually a random sample from some probability distribution, and we're trying to work out the parameters of that distribution. In greta we do that with the `distribution()` assignment function: ```{r likelihood} distribution(y) <- normal(mean, sd) ``` With the syntax `distribution(
### Defining the model Now all of the greta arrays making up the model have been created, we need to combine them and set up the model so that we can sample from it, using `model()`: ```{r hidden_model, echo = FALSE} x <- as_data(iris$Petal.Length) y <- as_data(iris$Sepal.Length) int <- normal(0, 1) coef <- normal(0, 3) sd <- student(3, 0, 1, truncation = c(0, Inf)) mean <- int + coef * x distribution(y) <- normal(mean, sd) ``` ```{r define_model} m <- model(int, coef, sd) ``` `model()` returns a 'greta model' object, which combines all of the greta arrays that make up the model. We can pass greta arrays as arguments to `model()` to flag them as the parameters we're interested in. When sampling from the model with `mcmc()` those will be the greta arrays for which samples will be returned. Alternatively, we can run `model()` without passing any greta arrays, in which case all of the greta arrays (except for data greta arrays) in the working environment will be set as the targets for sampling instead.
### Plotting greta provides a plot function for greta models to help you visualise and check the model before sampling from it. ```{r plot, eval = FALSE} plot(m) ``` ```{r plot_hidden, echo = FALSE} gr <- plot(m) DiagrammeR::export_graph(attr(gr, "dgr_graph"), file_name = "figures/full_graph.png", file_type = "png", width = 958 * 2, height = 450 * 2) ```
### Sampling When defining the model, greta combines all of the distributions together to define the *joint density* of the model, a measure of how 'good' (or how probable if we're being Bayesian) are a particular candidate set of values for the variables in the model. Now we have a greta model that will give us the joint density for a candidate set of values, so we can use that to carry out inference on the model. We do that using an Markov chain Monte Carlo (MCMC) algorithm to sample values of the parameters we're interested in, using the `mcmc()` function: ```{r mcmc, message = FALSE, results = "hide", progress = FALSE} draws <- mcmc(m, n_samples = 1000) ``` Here we're using 1000 steps of the static Hamiltonian Monte Carlo (HMC) algorithm on each of 4 separate chains, giving us 4000 samples. HMC uses the gradients of the joint density to efficiently explore the set of parameters. By default, greta also spends 1000 iterations on each chain 'warming up' (tuning the sampler parameters) and 'burning in' (moving to the area of highest probability) the sampler. `draws` is a `greta_mcmc_list` object, which inherits from the coda R package's `mcmc.list` object. So we can use functions from coda, or one of the many other MCMC software packages that use this format, to plot and summarise the MCMC samples. ```{r coda_summary} summary(draws) ``` The bayesplot package makes some nice plots of the MCMC chain and the parameter estimates ```{r mcmcvis, echo = TRUE, message = FALSE, out.width = c('400px', '400px'), fig.height=4, fig.width=5, fig.show='hold'} library (bayesplot) # set theme to avoid issues with fonts ggplot2::theme_set(ggplot2::theme_bw()) mcmc_trace(draws, facet_args = list(nrow = 3, ncol = 1)) mcmc_intervals(draws) ``` If your model is taking a long time whilst in the sampling phase and you want to take a look at the samples. You can stop the sampler (e.g. using the stop button in RStudio) and then retrieve the samples drawn so far, by using `stashed_samples()`. Note that this won't return anything if you stop the model during the warmup phase (since those aren't valid posterior samples) or if the sampler completed successfully. #### Tweaking the sampler greta's default sampler is (static) Hamiltonian Monte Carlo. The sampler will automatically tune itself during the warmup phase, to make it as efficient as possible. If the chain looks like it's moving too slowly, or if you are getting a lot of messages about proposals being rejected, the first thing to try is increasing the length of the warmup period from its default of 1000 iterations (via the `warmup` argument). If you're still getting a lot of rejected samples, it's often a good idea to manually set the initial values for the sampler (via the `initial_values` argument). This is often the case when you have lots of data; the more information there is, the more extreme the log probability, and the higher the risk of numerical problems. A downside of HMC is that it can't be used to sample discrete variables (e.g. integers), so we can't specify a model with a discrete distribution (e.g. Poisson or Binomial), unless it's in the likelihood. If we try to build such a model, greta will give us an error when we run `model()`. Future versions of greta will implement a wider range of MCMC samplers, including some for discrete variables.