diff --git a/DESCRIPTION b/DESCRIPTION
index a0df548c..f26486fd 100644
--- a/DESCRIPTION
+++ b/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 = "rich.fitzjohn@gmail.com"),
              person("Wes", "Hinsley", role = "aut"),
diff --git a/NAMESPACE b/NAMESPACE
index da753716..48ded14a 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/combine.R b/R/combine.R
index c25ce06b..de95e59d 100644
--- a/R/combine.R
+++ b/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()
diff --git a/R/dsl-differentiate-expr.R b/R/dsl-differentiate-expr.R
index fd1bdbc1..689a5a9e 100644
--- a/R/dsl-differentiate-expr.R
+++ b/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.
diff --git a/R/dsl-distributions.R b/R/dsl-distributions.R
index 3b7d8cfc..d8af6aae 100644
--- a/R/dsl-distributions.R
+++ b/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
 }
diff --git a/R/dsl-parse-error.R b/R/dsl-parse-error.R
index b9019376..3c08d7c8 100644
--- a/R/dsl-parse-error.R
+++ b/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)
 }
diff --git a/R/dsl.R b/R/dsl.R
index 4a6852b9..419bd224 100644
--- a/R/dsl.R
+++ b/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
diff --git a/R/example.R b/R/example.R
new file mode 100644
index 00000000..e15e46e9
--- /dev/null
+++ b/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))))
+}
diff --git a/R/model-function.R b/R/model-function.R
index 703ff4ba..af52af7e 100644
--- a/R/model-function.R
+++ b/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")
diff --git a/R/model.R b/R/model.R
index e10185e4..acd9663e 100644
--- a/R/model.R
+++ b/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])}")
+  }
+}
diff --git a/R/packer.R b/R/packer.R
index 32d9b056..7095c6b1 100644
--- a/R/packer.R
+++ b/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()
diff --git a/R/runner-simultaneous.R b/R/runner-simultaneous.R
index cdb5f88e..9a8f32b9 100644
--- a/R/runner-simultaneous.R
+++ b/R/runner-simultaneous.R
@@ -14,6 +14,11 @@
 ##'   [monty_sample()]
 ##'
 ##' @export
+##' @examples
+##' m <- monty_example("banana")
+##' s <- monty_sampler_random_walk(vcv = diag(2) * 0.01)
+##' r <- monty_runner_simultaneous()
+##' samples <- monty_sample(m, s, 200, runner = r)
 monty_runner_simultaneous <- function(progress = NULL) {
   validate_suitable <- function(model, observer) {
     require_multiple_parameters(
diff --git a/R/sample.R b/R/sample.R
index 32f468af..e13265e7 100644
--- a/R/sample.R
+++ b/R/sample.R
@@ -40,8 +40,12 @@
 ##'   restartable.  This will add additional data to the chains
 ##'   object.
 ##'
-##' @return A list of parameters and densities; we'll write tools for
-##'   dealing with this later.  Elements include:
+##' @return A list of parameters and densities.  We provide conversion
+##'   to formats used by other packages, notably
+##'   [posterior::as_draws_array], [posterior::as_draws_df] and
+##'   [coda::as.mcmc.list]; please let us know if you need conversion
+##'   to something else.  If you want to work directly with the
+##'   output, the elements in the list include:
 ##'
 ##' * `pars`: An array with three dimensions representing (in turn)
 ##'   parameter, sample and chain, so that `pars[i, j, k]` is the
@@ -65,6 +69,23 @@
 ##'   one is also subject to change.
 ##'
 ##' @export
+##' @examples
+##' m <- monty_example("banana")
+##' s <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
+##' samples <- monty_sample(m, s, 2000)
+##'
+##' # Quick conversion of parameters into something plottable:
+##' pars <- t(drop(samples$pars))
+##' plot(pars, pch = 19, cex = 0.75, col = "#0000ff55")
+##'
+##' # If you have the posterior package you might prefer converting to
+##' # its format for performing diagnoses:
+##' @examplesIf requireNamespace("posterior")
+##' res <- posterior::as_draws_df(samples)
+##' posterior::summarise_draws(res)
+##'
+##' # At this point you could also use the 'bayesplot' package to plot
+##' # diagnostics.
 monty_sample <- function(model, sampler, n_steps, initial = NULL,
                          n_chains = 1L, runner = NULL, observer = NULL,
                          restartable = FALSE) {
diff --git a/R/sampler-adaptive.R b/R/sampler-adaptive.R
index 514ef6c5..cf6c218c 100644
--- a/R/sampler-adaptive.R
+++ b/R/sampler-adaptive.R
@@ -96,6 +96,26 @@
 ##'   [monty_sample]
 ##'
 ##' @export
+##' @examples
+##' m <- monty_example("gaussian", matrix(c(1, 0.5, 0.5, 2), 2, 2))
+##' vcv <- diag(2) * 0.1
+##'
+##' # Sampling with a random walk
+##' s_rw <- monty_sampler_random_walk(vcv)
+##' res_rw <- monty_sample(m, s_rw, 1000)
+##'
+##' s_adapt <- monty_sampler_adaptive(vcv)
+##' res_adapt <- monty_sample(m, s_adapt, 1000)
+##'
+##' plot(drop(res_adapt$density), type = "l", col = 4)
+##' lines(drop(res_rw$density), type = "l", col = 2)
+##'
+##' # Estimated vcv from the sampler at the end of the simulation
+##' s_adapt$details[[1]]$vcv
+##'
+##' @examplesIf requireNamespace("coda")
+##' coda::effectiveSize(coda::as.mcmc.list(res_rw))
+##' coda::effectiveSize(coda::as.mcmc.list(res_adapt))
 monty_sampler_adaptive <- function(initial_vcv,
                                    initial_vcv_weight = 1000,
                                    initial_scaling = 1,
diff --git a/R/samples.R b/R/samples.R
index 8e9f1795..e28f1992 100644
--- a/R/samples.R
+++ b/R/samples.R
@@ -14,8 +14,11 @@ print.monty_samples <- function(x, ...) {
                   c("coda", "as.mcmc.list"))
   status <- suggested_package_status(target[, 1])
   status_cls <- c(missing = "danger", installed = "warning", loaded = "success")
-  target_str <- sprintf("{.alert-%s %s::%s() {cli::col_grey('[package %s]')}}",
-                        status_cls[status], target[, 1], target[, 2], status)
+  status_str <- status
+  status_str[status == "installed"] <- "installed, but not loaded"
+  target_str <- sprintf(
+    "{.alert-%s %s::%s() {cli::col_grey('[package %s]')}}",
+    status_cls[status], target[, 1], target[, 2], status_str)
   cli::cli_alert_info(
     "Conversion to other types is possible:")
   cli::cli_bullets(set_names(target_str, ">"))
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 3ed7b030..23dd80a0 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -17,6 +17,10 @@ reference:
     - monty_model_gradient
     - monty_model_direct_sample
 
+  - subtitle: Example
+    contents:
+    - monty_example
+
   - title: Domain specific language
     contents:
       - monty_dsl
diff --git a/man/monty_differentiation.Rd b/man/monty_differentiation.Rd
index 5008bc40..ec8e4290 100644
--- a/man/monty_differentiation.Rd
+++ b/man/monty_differentiation.Rd
@@ -50,3 +50,15 @@ We may need to make this slightly extensible in future, but for
 now the set of functions that can be differentiated is closed.
 }
 
+\section{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.
+}
+
+\examples{
+d <- monty_differentiation()
+d$differentiate(quote(sqrt(sin(x))), "x")
+D(quote(sqrt(sin(x))), "x")
+}
diff --git a/man/monty_dsl.Rd b/man/monty_dsl.Rd
index 57bc183f..352799fe 100644
--- a/man/monty_dsl.Rd
+++ b/man/monty_dsl.Rd
@@ -37,8 +37,7 @@ change name in future, as will its interface.
 }
 \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)
diff --git a/man/monty_dsl_distributions.Rd b/man/monty_dsl_distributions.Rd
index 64da993a..7a04b096 100644
--- a/man/monty_dsl_distributions.Rd
+++ b/man/monty_dsl_distributions.Rd
@@ -28,3 +28,6 @@ is primarily intended for use in packages which use
 information about which distributions and arguments would succeed
 there.
 }
+\examples{
+monty_dsl_distributions()
+}
diff --git a/man/monty_dsl_error_explain.Rd b/man/monty_dsl_error_explain.Rd
index 24727f2a..b18ad047 100644
--- a/man/monty_dsl_error_explain.Rd
+++ b/man/monty_dsl_error_explain.Rd
@@ -26,3 +26,6 @@ current implementation of this will send you to the rendered
 vignettes, but in the future we will arrange for offline rendering
 too.
 }
+\examples{
+monty_dsl_error_explain("E201")
+}
diff --git a/man/monty_dsl_parse_distribution.Rd b/man/monty_dsl_parse_distribution.Rd
index 9b020aa5..ebd80470 100644
--- a/man/monty_dsl_parse_distribution.Rd
+++ b/man/monty_dsl_parse_distribution.Rd
@@ -46,3 +46,10 @@ links through to the C++ API.  This function is designed for use
 from other packages that use monty, and is unlikely to be
 useful to most users.
 }
+\examples{
+# A successful match
+monty_dsl_parse_distribution(quote(Normal(0, 1)))
+
+# An unsuccessful match
+monty_dsl_parse_distribution(quote(Normal()))
+}
diff --git a/man/monty_example.Rd b/man/monty_example.Rd
new file mode 100644
index 00000000..f831f1c4
--- /dev/null
+++ b/man/monty_example.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/example.R
+\name{monty_example}
+\alias{monty_example}
+\title{Example models}
+\usage{
+monty_example(name, ...)
+}
+\arguments{
+\item{name}{Name of the example, as a string.  See Details for
+supported models.}
+
+\item{...}{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.}
+}
+\value{
+A \link{monty_model} object
+}
+\description{
+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.
+}
+\section{Supported models}{
+\subsection{\code{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 \code{alpha} and \code{beta} and is based on two
+successive draws, one conditional on the other.
+
+You can vary \code{sigma} for this model on creation, the default is 0.5
+}
+
+\subsection{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.
+}
+}
+
+\examples{
+monty_example("banana")
+monty_example("gaussian", diag(2))
+}
diff --git a/man/monty_model.Rd b/man/monty_model.Rd
index 5a89590e..a553a0cd 100644
--- a/man/monty_model.Rd
+++ b/man/monty_model.Rd
@@ -112,3 +112,7 @@ independent (for example, changing the parameters in group 2 does
 not affect the density of parameters proposed in group 3).
 }
 }
+\seealso{
+\link{monty_model_function}, which provides a simple interface
+for creating models from functions.
+}
diff --git a/man/monty_model_combine.Rd b/man/monty_model_combine.Rd
index 1fe628e7..1c8016c5 100644
--- a/man/monty_model_combine.Rd
+++ b/man/monty_model_combine.Rd
@@ -70,3 +70,21 @@ two elements corresponding to the first and second model (not the
 only part that makes a distinction between the two models here;
 for all components above they are equivalent.
 }
+\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))
+}
diff --git a/man/monty_model_density.Rd b/man/monty_model_density.Rd
index 88cdd9b3..ae606213 100644
--- a/man/monty_model_density.Rd
+++ b/man/monty_model_density.Rd
@@ -18,6 +18,11 @@ A log-density value, or vector of log-density values
 Compute log density for a model.  This is a wrapper around the
 \verb{$density} property within a \link{monty_model} object.
 }
+\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))
+}
 \seealso{
 \link{monty_model_gradient} for computing gradients and
 \link{monty_model_direct_sample} for sampling from a model.
diff --git a/man/monty_model_direct_sample.Rd b/man/monty_model_direct_sample.Rd
index 52b3992a..5eda0f4c 100644
--- a/man/monty_model_direct_sample.Rd
+++ b/man/monty_model_direct_sample.Rd
@@ -23,3 +23,14 @@ A vector or matrix of sampled parameters
 Directly sample from a model.  Not all models support this, and an
 error will be thrown if it is not possible.
 }
+\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)
+}
diff --git a/man/monty_model_function.Rd b/man/monty_model_function.Rd
index 9ddd5de9..7978b125 100644
--- a/man/monty_model_function.Rd
+++ b/man/monty_model_function.Rd
@@ -36,3 +36,18 @@ This interface will expand in future versions of monty to support
 gradients, stochastic models, parameter groups and simultaneous
 calculation of density.
 }
+\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))
+}
diff --git a/man/monty_model_gradient.Rd b/man/monty_model_gradient.Rd
index ecea6964..143d6c1d 100644
--- a/man/monty_model_gradient.Rd
+++ b/man/monty_model_gradient.Rd
@@ -22,6 +22,15 @@ Compute the gradient of log density (which is returned by
 \link{monty_model_density}) with respect to parameters.  Not all models
 support this, and an error will be thrown if it is not possible.
 }
+\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))
+}
 \seealso{
 \link{monty_model_density} for log density, and
 \link{monty_model_direct_sample} to sample from a model
diff --git a/man/monty_model_properties.Rd b/man/monty_model_properties.Rd
index 5e67a945..a7a02b90 100644
--- a/man/monty_model_properties.Rd
+++ b/man/monty_model_properties.Rd
@@ -54,3 +54,10 @@ but you can pass the return value of this as the \code{properties}
 argument of monty_model to enforce that your model does actually
 have these properties.
 }
+\examples{
+# Default properties:
+monty_model_properties()
+
+# Set some properties:
+monty_model_properties(has_gradient = TRUE, is_stochastic = FALSE)
+}
diff --git a/man/monty_packer.Rd b/man/monty_packer.Rd
index 8affac6e..7ae2ff25 100644
--- a/man/monty_packer.Rd
+++ b/man/monty_packer.Rd
@@ -125,8 +125,8 @@ vector in order to use within an MCMC or optimisation algorithm.
 Ultimately, our unpacked vector will need to hold four elements
 (\code{b11}, \code{b12}, \code{b21}, \code{b22}), but there are only three distinct
 values as the two off-diagonal elements will be the same (i.e.,
-\verb{b12 == b21``).  So we might write this passing in }b_raw = 3\code{to}array\verb{, so that our unpacked list holds }b_raw = c(b11, b12,
-b22)\verb{.  We would then write }process` as something like:
+\code{b12 == b21}).  So we might write this passing in \code{b_raw = 3} to
+\code{array}, so that our unpacked list holds \code{b_raw = c(b11, b12, b22)}.  We would then write \code{process} as something like:
 
 \if{html}{\out{<div class="sourceCode">}}\preformatted{process <- function(x) \{
   list(b = matrix(x$b_raw[c(1, 2, 2, 3)], 2, 2))
@@ -172,3 +172,46 @@ back up, though it's not hard.  Please let us know if you would
 use this.
 }
 
+\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))
+}
diff --git a/man/monty_runner_simultaneous.Rd b/man/monty_runner_simultaneous.Rd
index e9bb3613..2eea2188 100644
--- a/man/monty_runner_simultaneous.Rd
+++ b/man/monty_runner_simultaneous.Rd
@@ -28,3 +28,9 @@ be faster than running in parallel, but primarily this exists so
 that we can see that samplers can work with multiple samples at
 once.
 }
+\examples{
+m <- monty_example("banana")
+s <- monty_sampler_random_walk(vcv = diag(2) * 0.01)
+r <- monty_runner_simultaneous()
+samples <- monty_sample(m, s, 200, runner = r)
+}
diff --git a/man/monty_sample.Rd b/man/monty_sample.Rd
index 546a53e2..a7c1ef5a 100644
--- a/man/monty_sample.Rd
+++ b/man/monty_sample.Rd
@@ -50,8 +50,12 @@ restartable.  This will add additional data to the chains
 object.}
 }
 \value{
-A list of parameters and densities; we'll write tools for
-dealing with this later.  Elements include:
+A list of parameters and densities.  We provide conversion
+to formats used by other packages, notably
+\link[posterior:draws_array]{posterior::as_draws_array}, \link[posterior:draws_df]{posterior::as_draws_df} and
+\link[coda:mcmc.list]{coda::as.mcmc.list}; please let us know if you need conversion
+to something else.  If you want to work directly with the
+output, the elements in the list include:
 \itemize{
 \item \code{pars}: An array with three dimensions representing (in turn)
 parameter, sample and chain, so that \code{pars[i, j, k]} is the
@@ -79,3 +83,22 @@ support for distributing over workers, and for things like
 parallel reproducible streams of random numbers.  For now it just
 runs a single chain as a proof of concept.
 }
+\examples{
+m <- monty_example("banana")
+s <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
+samples <- monty_sample(m, s, 2000)
+
+# Quick conversion of parameters into something plottable:
+pars <- t(drop(samples$pars))
+plot(pars, pch = 19, cex = 0.75, col = "#0000ff55")
+
+# If you have the posterior package you might prefer converting to
+# its format for performing diagnoses:
+\dontshow{if (requireNamespace("posterior")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+res <- posterior::as_draws_df(samples)
+posterior::summarise_draws(res)
+
+# At this point you could also use the 'bayesplot' package to plot
+# diagnostics.
+\dontshow{\}) # examplesIf}
+}
diff --git a/man/monty_sampler_adaptive.Rd b/man/monty_sampler_adaptive.Rd
index 79014dcd..3656d8a4 100644
--- a/man/monty_sampler_adaptive.Rd
+++ b/man/monty_sampler_adaptive.Rd
@@ -117,3 +117,25 @@ Spencer SEF (2021) Accelerating adaptation in the adaptive
 Metropolis–Hastings random walk algorithm. Australian & New Zealand Journal
 of Statistics 63:468-484.
 }
+\examples{
+m <- monty_example("gaussian", matrix(c(1, 0.5, 0.5, 2), 2, 2))
+vcv <- diag(2) * 0.1
+
+# Sampling with a random walk
+s_rw <- monty_sampler_random_walk(vcv)
+res_rw <- monty_sample(m, s_rw, 1000)
+
+s_adapt <- monty_sampler_adaptive(vcv)
+res_adapt <- monty_sample(m, s_adapt, 1000)
+
+plot(drop(res_adapt$density), type = "l", col = 4)
+lines(drop(res_rw$density), type = "l", col = 2)
+
+# Estimated vcv from the sampler at the end of the simulation
+s_adapt$details[[1]]$vcv
+
+\dontshow{if (requireNamespace("coda")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+coda::effectiveSize(coda::as.mcmc.list(res_rw))
+coda::effectiveSize(coda::as.mcmc.list(res_adapt))
+\dontshow{\}) # examplesIf}
+}
diff --git a/tests/testthat/helper-monty.R b/tests/testthat/helper-monty.R
index 9acd5004..3c8efbcb 100644
--- a/tests/testthat/helper-monty.R
+++ b/tests/testthat/helper-monty.R
@@ -176,53 +176,6 @@ ex_dust_sir <- function(n_particles = 100, n_threads = 1,
 }
 
 
-ex_simple_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))))
-}
-
-
-ex_banana <- function(sd = 0.5) {
-  monty_model(
-    list(parameters = c("a", "b"),
-         direct_sample = function(rng) {
-           b <- rng$random_normal(1)
-           a <- rng$normal(1, b^2, sd)
-           c(a, b)
-         },
-         density = function(x) {
-           if (length(dim2(x)) == 1) {
-             a <- x[1]
-             b <- x[2]
-           } else {
-             a <- x[1, ]
-             b <- x[2, ]
-           }
-           dnorm(b, log = TRUE) + dnorm((a - b^2) / sd, log = TRUE)
-         },
-         gradient = function(x) {
-           if (length(dim2(x)) == 1) {
-             a <- x[1]
-             b <- x[2]
-             c((b^2 - a) / sd^2,
-               -b + 2 * b * (a - b^2) / sd^2)
-           } else {
-             a <- x[1, ]
-             b <- x[2, ]
-             rbind((b^2 - a) / sd^2,
-                   -b + 2 * b * (a - b^2) / sd^2)
-           }
-         },
-         domain = cbind(rep(-Inf, 2), rep(Inf, 2))),
-    monty_model_properties(allow_multiple_parameters = TRUE))
-}
-
-
 random_array <- function(dim, named = FALSE) {
   if (named) {
     dn <- lapply(seq_along(dim), function(i) {
diff --git a/tests/testthat/test-example.R b/tests/testthat/test-example.R
new file mode 100644
index 00000000..e8b5dd96
--- /dev/null
+++ b/tests/testthat/test-example.R
@@ -0,0 +1,5 @@
+test_that("error if unexpected example used", {
+  expect_error(
+    monty_example("unknown"),
+    "'name' must be one of")
+})
diff --git a/tests/testthat/test-model-wrappers.R b/tests/testthat/test-model-wrappers.R
index b1b5880c..ee0005cf 100644
--- a/tests/testthat/test-model-wrappers.R
+++ b/tests/testthat/test-model-wrappers.R
@@ -73,14 +73,14 @@ test_that("can directly sample from a model", {
 
 
 test_that("can compute multiple gradients", {
-  m <- ex_banana()
+  m <- monty_example("banana")
   expect_equal(monty_model_gradient(m, c(0, 0)), c(0, 0))
   expect_equal(monty_model_gradient(m, c(0, 0), named = TRUE),
-               c(a = 0, b = 0))
+               c(alpha = 0, beta = 0))
 
   p <- cbind(c(0, 0), c(1, 0), c(0, 1))
   expect_equal(monty_model_gradient(m, p),
                rbind(c(0, -4, 4), c(0, 0, -9)))
   expect_equal(monty_model_gradient(m, p, named = TRUE),
-               rbind(a = c(0, -4, 4), b = c(0, 0, -9)))
+               rbind(alpha = c(0, -4, 4), beta = c(0, 0, -9)))
 })
diff --git a/tests/testthat/test-sampler-adaptive.R b/tests/testthat/test-sampler-adaptive.R
index 0d745561..e6d9d11d 100644
--- a/tests/testthat/test-sampler-adaptive.R
+++ b/tests/testthat/test-sampler-adaptive.R
@@ -1,5 +1,5 @@
 test_that("Empirical VCV calculated correctly with forget_rate = 0", {
-  m <- ex_simple_gaussian(vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
+  m <- monty_example("gaussian", vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
 
   sampler <- monty_sampler_adaptive(initial_vcv = diag(c(0.01, 0.01)),
                                     forget_rate = 0,
@@ -17,7 +17,7 @@ test_that("Empirical VCV calculated correctly with forget_rate = 0", {
 
 
 test_that("Empirical VCV calculated correctly with forget_rate = 0.1", {
-  m <- ex_simple_gaussian(vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
+  m <- monty_example("gaussian", vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
 
   sampler <- monty_sampler_adaptive(initial_vcv = diag(c(0.01, 0.01)),
                                     forget_rate = 0.1)
@@ -35,7 +35,7 @@ test_that("Empirical VCV calculated correctly with forget_rate = 0.1", {
 
 
 test_that("Empirical VCV correct using both forget_rate and forget_end", {
-  m <- ex_simple_gaussian(vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
+  m <- monty_example("gaussian", vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
 
   sampler <- monty_sampler_adaptive(initial_vcv = diag(c(0.01, 0.01)),
                                     forget_rate = 0.5,
@@ -55,7 +55,7 @@ test_that("Empirical VCV correct using both forget_rate and forget_end", {
 
 
 test_that("Empirical VCV correct using forget_rate, forget_end and adapt_end", {
-  m <- ex_simple_gaussian(vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
+  m <- monty_example("gaussian", vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
 
   sampler <- monty_sampler_adaptive(initial_vcv = diag(c(0.01, 0.01)),
                                     forget_rate = 0.25,
@@ -76,7 +76,7 @@ test_that("Empirical VCV correct using forget_rate, forget_end and adapt_end", {
 
 
 test_that("can continue adaptive sampler", {
-  m <- ex_simple_gaussian(vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
+  m <- monty_example("gaussian", vcv = rbind(c(0.02, 0.01), c(0.01, 0.03)))
   sampler <- monty_sampler_adaptive(initial_vcv = diag(c(0.01, 0.01)))
 
   set.seed(1)
diff --git a/tests/testthat/test-sampler-hmc.R b/tests/testthat/test-sampler-hmc.R
index 8b153f20..b39b13f2 100644
--- a/tests/testthat/test-sampler-hmc.R
+++ b/tests/testthat/test-sampler-hmc.R
@@ -1,5 +1,5 @@
 test_that("can run hmc", {
-  m <- ex_simple_gaussian(diag(2))
+  m <- monty_example("gaussian", diag(2))
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
   set.seed(1)
   res <- monty_sample(m, sampler, 2000)
@@ -10,7 +10,7 @@ test_that("can run hmc", {
 
 
 test_that("can run hmc on banana shaped model", {
-  m <- ex_banana()
+  m <- monty_example("banana")
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
   set.seed(1)
   res <- monty_sample(m, sampler, 2000)
@@ -23,7 +23,7 @@ test_that("can run hmc on banana shaped model", {
 
 
 test_that("can output debug traces", {
-  m <- ex_simple_gaussian(diag(2))
+  m <- monty_example("gaussian", diag(2))
 
   set.seed(1)
   sampler1 <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
@@ -49,7 +49,7 @@ test_that("can output debug traces", {
 
 
 test_that("can use custom vcv", {
-  m <- ex_simple_gaussian(diag(2))
+  m <- monty_example("gaussian", diag(2))
   set.seed(1)
   sampler1 <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
   res1 <- monty_sample(m, sampler1, 30)
@@ -69,7 +69,7 @@ test_that("can use custom vcv", {
 
 
 test_that("given vcv and parameters must be compatible", {
-  m <- ex_simple_gaussian(diag(2))
+  m <- monty_example("gaussian", diag(2))
   set.seed(1)
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10,
                                  vcv = diag(3))
@@ -174,7 +174,7 @@ test_that("prevent weird distributions", {
 
 
 test_that("can continue hmc without debug", {
-  m <- ex_simple_gaussian(diag(2))
+  m <- monty_example("gaussian", diag(2))
 
   set.seed(1)
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
@@ -189,7 +189,7 @@ test_that("can continue hmc without debug", {
 
 
 test_that("can continue hmc with debug", {
-  m <- ex_simple_gaussian(diag(2))
+  m <- monty_example("gaussian", diag(2))
 
   set.seed(1)
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10,
@@ -225,7 +225,7 @@ test_that("can't use hmc with stochastic models", {
 
 
 test_that("can run hmc model simultaneously", {
-  m <- ex_banana()
+  m <- monty_example("banana")
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10)
   runner <- monty_runner_simultaneous()
   set.seed(1)
@@ -237,7 +237,7 @@ test_that("can run hmc model simultaneously", {
 
 
 test_that("can run hmc model simultaneously, with debug", {
-  m <- ex_banana()
+  m <- monty_example("banana")
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10,
                                  debug = TRUE)
   runner <- monty_runner_simultaneous()
@@ -253,7 +253,7 @@ test_that("can run hmc model simultaneously, with debug", {
 
 
 test_that("can continue a hmc model simultaneously, with debug", {
-  m <- ex_banana()
+  m <- monty_example("banana")
   sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10,
                                  debug = TRUE)
   runner <- monty_runner_simultaneous()
diff --git a/tests/testthat/test-sampler-nested-adaptive.R b/tests/testthat/test-sampler-nested-adaptive.R
index 439f0beb..cd586a56 100644
--- a/tests/testthat/test-sampler-nested-adaptive.R
+++ b/tests/testthat/test-sampler-nested-adaptive.R
@@ -46,7 +46,7 @@ test_that("validate vcv inputs on construction of sampler", {
 
 
 test_that("can't use nested sampler with models that are not nested", {
-  m <- ex_simple_gaussian(diag(3))
+  m <- monty_example("gaussian", diag(3))
   vcv <- list(base = diag(1), groups = list(diag(2), diag(3)))
   s <- monty_sampler_nested_adaptive(vcv)
   expect_error(
diff --git a/tests/testthat/test-sampler-nested-random-walk.R b/tests/testthat/test-sampler-nested-random-walk.R
index 8b22fa9c..7d070250 100644
--- a/tests/testthat/test-sampler-nested-random-walk.R
+++ b/tests/testthat/test-sampler-nested-random-walk.R
@@ -46,7 +46,7 @@ test_that("validate vcv inputs on construction of sampler", {
 
 
 test_that("can't use nested sampler with models that are not nested", {
-  m <- ex_simple_gaussian(diag(3))
+  m <- monty_example("gaussian", diag(3))
   vcv <- list(base = diag(1), groups = list(diag(2), diag(3)))
   s <- monty_sampler_nested_random_walk(vcv)
   expect_error(
diff --git a/vignettes/samplers.Rmd b/vignettes/samplers.Rmd
index 5d5f1b5a..6638cc21 100644
--- a/vignettes/samplers.Rmd
+++ b/vignettes/samplers.Rmd
@@ -37,24 +37,24 @@ with $\sigma$ the standard deviation of the conditional draw.
 We'll use R's `dnorm` for the likelihood calculations and then differentiate the density by hand (using the formula of the gaussian density) to derive the gradient by hand.  We skip the details here as this will be automated in an upcoming version of the package.
 
 ```{r}
-banana_model <- function(sd = 0.5) {
+banana_model <- 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, sd)
+      alpha <- rng$normal(1, beta^2, sigma)
       c(alpha, beta)
     },
     density = function(x) {
       alpha <- x[[1]]
       beta <- x[[2]]
-      dnorm(beta, log = TRUE) + dnorm((alpha - beta^2) / sd, log = TRUE)
+      dnorm(beta, log = TRUE) + dnorm((alpha - beta^2) / sigma, log = TRUE)
     },
     gradient = function(x) {
       alpha <- x[[1]]
       beta <- x[[2]]
-      c((beta^2 - alpha) / sd^2,
-        -beta + 2 * beta * (alpha - beta^2) / sd^2)
+      c((beta^2 - alpha) / sigma^2,
+        -beta + 2 * beta * (alpha - beta^2) / sigma^2)
     },
     domain = rbind(c(-Inf, Inf), c(-Inf, Inf))))
 }
@@ -89,11 +89,11 @@ It is also possible to compute the 95% confidence interval of the distribution u
 ```{r}
 theta <- seq(0, 2 * pi, length.out = 10000)
 z95 <- local({
-  sd <- 0.5
+  sigma <- 0.5
   r <- sqrt(qchisq(.95, df = 2))
   x <- r * cos(theta)
   y <- r * sin(theta)
-  cbind(x^2 + y * sd, x)
+  cbind(x^2 + y * sigma, x)
 })
 image(a, b, z, xlab = "alpha", ylab = "beta")
 lines(z95[, 1], z95[, 2])