Merge pull request #70 from mrc-ide/mrc-5765

Add basic examples

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: monty
 Title: Monte Carlo Models
-Version: 0.2.4
+Version: 0.2.5
 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
                     email = "[email protected]"),
              person("Wes", "Hinsley", role = "aut"),

NAMESPACE

@@ -6,6 +6,7 @@ S3method(coda::as.mcmc.list,monty_samples)
 S3method(posterior::as_draws_array,monty_samples)
 S3method(posterior::as_draws_df,monty_samples)
 S3method(print,monty_model)
+S3method(print,monty_model_properties)
 S3method(print,monty_observer)
 S3method(print,monty_packer)
 S3method(print,monty_runner)
@@ -17,6 +18,7 @@ export(monty_dsl)
 export(monty_dsl_distributions)
 export(monty_dsl_error_explain)
 export(monty_dsl_parse_distribution)
+export(monty_example)
 export(monty_model)
 export(monty_model_combine)
 export(monty_model_density)

R/combine.R

@@ -65,6 +65,23 @@
 ##' @return A [monty_model] object
 ##'
 ##' @export
+##' @examples
+##' # A simple example; a model that contains something of interest,
+##' # and a simple prior from monty_dsl
+##' likelihood <- monty_example("banana")
+##' prior <- monty_dsl({
+##'   alpha ~ Normal(0, 1)
+##'   beta ~ Normal(0, 10)
+##' })
+##' posterior <- likelihood + prior
+##' posterior
+##'
+##' # The same thing, more explicitly:
+##' monty_model_combine(likelihood, prior)
+##'
+##' # Control properties of the combined model:
+##' monty_model_combine(likelihood, prior,
+##'                     monty_model_properties(has_gradient = FALSE))
 monty_model_combine <- function(a, b, properties = NULL,
                                 name_a = "a", name_b = "b") {
   call <- environment()

R/dsl-differentiate-expr.R

@@ -26,6 +26,13 @@
 ##' We may need to make this slightly extensible in future, but for
 ##' now the set of functions that can be differentiated is closed.
 ##'
+##' # Warning
+##'
+##' The way of accessing distribution support here is peculiar and the
+##' return type unusual.  This is intentional, and we expect a more
+##' conventional interface in the future once this package settles
+##' down.
+##'
 ##' @title Differentiate expressions
 ##'
 ##' @return A list of related objects:
@@ -42,6 +49,10 @@
 ##' to differentiate to allow programs to fail fast.
 ##'
 ##' @export
+##' @examples
+##' d <- monty_differentiation()
+##' d$differentiate(quote(sqrt(sin(x))), "x")
+##' D(quote(sqrt(sin(x))), "x")
 monty_differentiation <- function() {
   ## TODO: we might want to export a stripped down copy of maths
   ## perhaps.

R/dsl-distributions.R

@@ -20,6 +20,8 @@
 ##'   supporting everything this way soon.
 ##'
 ##' @export
+##' @examples
+##' monty_dsl_distributions()
 monty_dsl_distributions <- function() {
   dsl_distribution_summary
 }

R/dsl-parse-error.R

@@ -18,6 +18,8 @@
 ##' @return Nothing, this is called for its side effect only
 ##'
 ##' @export
+##' @examples
+##' monty_dsl_error_explain("E201")
 monty_dsl_error_explain <- function(code, how = "pretty") {
   error_explain(dsl_errors, code, how)
 }

R/dsl.R

@@ -29,8 +29,7 @@
 ##' @export
 ##' @examples
 ##'
-##' # Expressions that correspond to models can be passed in with no
-##' # quoting
+##' # Expressions that define models can be passed in with no quoting
 ##' monty_dsl(a ~ Normal(0, 1))
 ##' monty_dsl({
 ##'   a ~ Normal(0, 1)
@@ -109,6 +108,12 @@ monty_dsl_parse <- function(x, type = NULL, gradient = NULL) {
 ##'   a first argument) a rng object (see [monty_rng])
 ##'
 ##' @export
+##' @examples
+##' # A successful match
+##' monty_dsl_parse_distribution(quote(Normal(0, 1)))
+##'
+##' # An unsuccessful match
+##' monty_dsl_parse_distribution(quote(Normal()))
 monty_dsl_parse_distribution <- function(expr, name = NULL) {
   ## Here, the user has not provided a call to anything, or a call to
   ## something that is not recognised as a distribution.  We throw the

R/example.R

@@ -0,0 +1,99 @@
+##' Load example models from monty.  These models exist so that we can
+##' create (hopefully) interesting examples in the documentation
+##' without them becoming overwhelming.  You should probably not use
+##' these for anything other than exploring the package.
+##'
+##' # Supported models
+##'
+##' ## `bananna`
+##'
+##' The banana model is a two-dimensional banana-shaped function,
+##' picked because it is quite annoying to sample from directly.  The
+##' model has two parameters `alpha` and `beta` and is based on two
+##' successive draws, one conditional on the other.
+##'
+##' You can vary `sigma` for this model on creation, the default is 0.5
+##'
+##' ## Gaussian
+##'
+##' A multivariate Gaussian centred at the origin.  Takes a
+##' variance-covariance-matrix as its argument.  Parameters are
+##' letters a, b, ... up to the number of dimensions.
+##'
+##' @title Example models
+##'
+##' @param name Name of the example, as a string.  See Details for
+##'   supported models.
+##'
+##' @param ... Optional parameters that are passed to create the
+##'   model.  All models can be created with no additional parameters,
+##'   but you can tweak their behaviour by passing named parameters
+##'   here.  See Details.
+##'
+##' @return A [monty_model] object
+##' @export
+##'
+##' @examples
+##' monty_example("banana")
+##' monty_example("gaussian", diag(2))
+monty_example <- function(name, ...) {
+  examples <- list(banana = monty_example_banana,
+                   gaussian = monty_example_gaussian)
+  examples[[match_value(name, names(examples))]](...)
+}
+
+
+monty_example_banana <- function(sigma = 0.5) {
+  monty_model(
+    list(
+      parameters = c("alpha", "beta"),
+      direct_sample = function(rng) {
+        beta <- rng$random_normal(1)
+        alpha <- rng$normal(1, beta^2, sigma)
+        if (length(alpha) == 1) {
+          c(alpha, beta)
+        } else {
+          rbind(alpha, beta)
+        }
+      },
+      density = function(x) {
+        if (is.matrix(x)) {
+          alpha <- x[1, ]
+          beta <- x[2, ]
+        } else {
+          alpha <- x[[1]]
+          beta <- x[[2]]
+        }
+        stats::dnorm(beta, log = TRUE) +
+          stats::dnorm((alpha - beta^2) / sigma, log = TRUE)
+      },
+      gradient = function(x) {
+        if (is.matrix(x)) {
+          alpha <- x[1, ]
+          beta <- x[2, ]
+        } else {
+          alpha <- x[[1]]
+          beta <- x[[2]]
+        }
+        da <- (beta^2 - alpha) / sigma^2
+        db <- -beta + 2 * beta * (alpha - beta^2) / sigma^2
+        if (is.matrix(x)) {
+          rbind(da, db, deparse.level = 0)
+        } else {
+          c(da, db)
+        }
+      },
+      domain = rbind(c(-Inf, Inf), c(-Inf, Inf))),
+    properties = monty_model_properties(allow_multiple_parameters = TRUE))
+}
+
+
+monty_example_gaussian <- function(vcv) {
+  n <- nrow(vcv)
+  monty_model(list(
+    parameters = letters[seq_len(n)],
+    direct_sample = make_rmvnorm(vcv, centred = TRUE),
+    density = make_ldmvnorm(vcv),
+    gradient = make_deriv_ldmvnorm(vcv),
+    domain = cbind(rep(-Inf, n), rep(Inf, n))))
+}

R/model-function.R

@@ -27,6 +27,20 @@
 ##'   representing all parameters.
 ##'
 ##' @export
+##' @examples
+##' banana <- function(a, b, sd) {
+##'   dnorm(b, log = TRUE) + dnorm((a - b^2) / sd, log = TRUE)
+##' }
+##' m <- monty_model_function(banana, fixed = list(sd = 0.25))
+##' m
+##'
+##' # Density from our new model. Note that this computes density
+##' # using an unstructured parameter vector, which is mapped to 'a'
+##' # and 'b':
+##' monty_model_density(m, c(0, 0))
+##'
+##' # Same as the built-in banana example:
+##' monty_model_density(monty_example("banana"), c(0, 0))
 monty_model_function <- function(density, packer = NULL, fixed = NULL) {
   if (!is.function(density)) {
     cli::cli_abort("Expected 'density' to be a function", arg = "density")

R/model.R

@@ -40,6 +40,13 @@
 ##'   not be modified.
 ##'
 ##' @export
+##'
+##' @examples
+##' # Default properties:
+##' monty_model_properties()
+##'
+##' # Set some properties:
+##' monty_model_properties(has_gradient = TRUE, is_stochastic = FALSE)
 monty_model_properties <- function(has_gradient = NULL,
                                    has_direct_sample = NULL,
                                    is_stochastic = NULL,
@@ -165,6 +172,9 @@ monty_model_properties <- function(has_gradient = NULL,
 ##'     * `has_parameter_groups`: The model has separable parameter groups
 ##'
 ##' @export
+##'
+##' @seealso [monty_model_function], which provides a simple interface
+##'   for creating models from functions.
 monty_model <- function(model, properties = NULL) {
   call <- environment() # for nicer stack traces
   parameters <- validate_model_parameters(model, call)
@@ -212,6 +222,10 @@ monty_model <- function(model, properties = NULL) {
 ##'   [monty_model_direct_sample] for sampling from a model.
 ##'
 ##' @export
+##' @examples
+##' m <- monty_model_function(function(a, b) dnorm(0, a, b, log = TRUE))
+##' monty_model_density(m, c(0, 1))
+##' monty_model_density(m, c(0, 10))
 monty_model_density <- function(model, parameters) {
   require_monty_model(model)
   check_model_parameters(model, parameters)
@@ -236,6 +250,14 @@ monty_model_density <- function(model, parameters) {
 ##'   [monty_model_direct_sample] to sample from a model
 ##'
 ##' @export
+##' @examples
+##' m <- monty_example("banana")
+##' # Global maximum at (0, 0), and gradient is zero there:
+##' monty_model_density(m, c(0, 0))
+##' monty_model_gradient(m, c(0, 0))
+##'
+##' # Nonzero gradient away from the origin:
+##' monty_model_gradient(m, c(0.4, 0.2))
 monty_model_gradient <- function(model, parameters, named = FALSE) {
   require_monty_model(model)
   require_gradient(
@@ -270,6 +292,16 @@ monty_model_gradient <- function(model, parameters, named = FALSE) {
 ##' @return A vector or matrix of sampled parameters
 ##'
 ##' @export
+##' @examples
+##' m <- monty_example("banana")
+##'
+##' r <- monty_rng$new()
+##' monty_model_direct_sample(m, r)
+##' monty_model_direct_sample(m, r, named = TRUE)
+##'
+##' r <- monty_rng$new(n_streams = 3)
+##' monty_model_direct_sample(m, r)
+##' monty_model_direct_sample(m, r, named = TRUE)
 monty_model_direct_sample <- function(model, rng, named = FALSE) {
   require_monty_model(model)
   require_direct_sample(
@@ -278,7 +310,11 @@ monty_model_direct_sample <- function(model, rng, named = FALSE) {
   assert_scalar_logical(named, call = environment())
   ret <- model$direct_sample(rng)
   if (named) {
-    names(ret) <- model$parameters
+    if (is.matrix(ret)) {
+      rownames(ret) <- model$parameters
+    } else {
+      names(ret) <- model$parameters
+    }
   }
   ret
 }
@@ -583,3 +619,20 @@ check_model_parameters <- function(model, parameters, call = parent.frame()) {
       arg = "parameters", call = call)
   }
 }
+
+
+##' @export
+print.monty_model_properties <- function(x, ...) {
+  cli::cli_h1("<monty_model_properties>")
+  unset <- vlapply(x, is.null)
+  is_set <- !unset
+  if (any(is_set)) {
+    cli::cli_bullets(
+      set_names(sprintf("%s: {.code %s}",
+                        names(x)[is_set], vcapply(x[is_set], as.character)),
+                "*"))
+  }
+  if (any(unset)) {
+    cli::cli_alert_info("Unset: {squote(names(x)[unset])}")
+  }
+}

R/packer.R

@@ -68,7 +68,7 @@
 ##' Ultimately, our unpacked vector will need to hold four elements
 ##' (`b11`, `b12`, `b21`, `b22`), but there are only three distinct
 ##' values as the two off-diagonal elements will be the same (i.e.,
-##' `b12 == b21``).  So we might write this passing in `b_raw = 3` to
+##' `b12 == b21`).  So we might write this passing in `b_raw = 3` to
 ##' `array`, so that our unpacked list holds `b_raw = c(b11, b12,
 ##' b22)`.  We would then write `process` as something like:
 ##'
@@ -168,6 +168,49 @@
 ##'   can be found.
 ##'
 ##' @export
+##'
+##' @examples
+##' # Here's a really simple example
+##' p <- monty_packer(c("a", "b", "c"))
+##' p
+##'
+##' p$pack(list(a = 1, b = 2, c = 3))
+##' p$unpack(1:3)
+##'
+##' # Sometimes we have a vector embedded in our parameters:
+##' p <- monty_packer(c("a", "b"), list(v = 4))
+##' p$pack(list(a = 1, b = 2, v = c(6, 7, 8, 9)))
+##' p$unpack(c(1, 2, 6, 7, 8, 9))
+##'
+##' # Or a higher dimensional structure such as a matrix:
+##' p <- monty_packer(c("a", "b"), list(m = c(2, 2)))
+##' p$unpack(c(1, 2, 6, 7, 8, 9))
+##'
+##' # You can use a packer to set "fixed" parameters that do not vary
+##' # with the underlying model being fit, but are required by your model.
+##' # This is simpler than the "closure" approach used previously in our
+##' # mcstate package and also easier to accommodate with differentiable
+##' # models:
+##' p <- monty_packer(
+##'   c("a", "b"),
+##'   fixed = list(d = data.frame(n = 1:3, m = runif(3))))
+##' p$unpack(1:2)
+##' p$pack(p$unpack(1:2))
+##'
+##' # The example from above, where we create a symmetric 2 x 2 matrix
+##' # from a 3-element vector, alongside a scalar:
+##' p <- monty_packer(
+##'   scalar = "a",
+##'   array = list(b_flat = 3),
+##'   process = function(p) list(b = matrix(p$b_flat[c(1, 2, 2, 3)], 2, 2)))
+##'
+##' # Unpacking we see "b_flat" is still in the list, but "b" is our
+##' # symmetric matrix:
+##' p$unpack(1:4)
+##'
+##' # The processed elements are ignored on the return pack:
+##' p$pack(list(a = 1, b_flat = 2:4, b = matrix(c(2, 3, 3, 4), 2, 2)))
+##' p$pack(list(a = 1, b_flat = 2:4))
 monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL,
                          process = NULL) {
   call <- environment()