Package 'monty'

Description

Differentiate expressions in the monty DSL. This function is exported for advanced use, and really so that we can use it from odin. But it has the potential to be generally useful, so while we'll tweak the interface quite a lot over the next while it is fine to use if you can handle some disruption.

Usage

monty_differentiation()

Details

R already has support for differentiating expressions using D, which is useful for creating derivatives of simple functions to pass into non-linear optimisation. We need something a bit more flexible for differentiating models in the monty DSL (monty_dsl) and also in the related odin DSL.

Value

A list of related objects:

• differentiate: A function that can differentiate an expression with respect to a variable (as a string).

• maths: Some mathematical utilities for constructing expressions. This will be documented later, but the most useful bits on here are the function elements times, plus and plus_fold.

We will expand this soon to advertise what functions we are able to differentiate to allow programs to fail fast.

Differences to D()

• We try a little harder to simplify expressions.

• The distribution functions in the monty DSL (e.g., Poisson) are (will be) handled specially, allowing substitution of log-densities and expectations.

• Once we support array expressions, we will be able to differentiate through these.

Roadmap

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.

Examples

d <- monty_differentiation()
d$differentiate(quote(sqrt(sin(x))), "x")
D(quote(sqrt(sin(x))), "x")

Description

Check and expand a domain, where it is used alongside a monty_packer object. This can be used to expand domains for logical parameters (e.g. a vector b) into its specific names (e.g., b[1], b[2], etc) without having to rely on the internals about how these names are constructed.

Usage

monty_domain_expand(domain, packer)

Arguments

Value

A two dimensional matrix representing your domain, or NULL if domain was given as NULL.

Examples

packer <- monty_packer(c("a", "b"), list(x = 3, y = c(2, 2)))
monty_domain_expand(NULL, packer)
monty_domain_expand(rbind(x = c(0, 1)), packer)
monty_domain_expand(rbind(x = c(0, 1), "x[2]" = c(0, Inf)), packer)
monty_domain_expand(rbind(x = c(0, 1), "y" = c(0, Inf)), packer)

Description

Create a model using the monty DSL; this function will likely change name in future, as will its interface.

Usage

monty_dsl(x, type = NULL, gradient = NULL, fixed = NULL)

Arguments

Value

A monty_model object derived from the expressions you provide.

Examples

# Expressions that define models can be passed in with no quoting
monty_dsl(a ~ Normal(0, 1))
monty_dsl({
  a ~ Normal(0, 1)
  b ~ Exponential(1)
})

# You can also pass strings
monty_dsl("a ~ Normal(0, 1)")

Description

Report information about supported distributions in the DSL. This is primarily intended for use in packages which use monty_dsl_parse_distribution, as this function reports information about which distributions and arguments would succeed there.

Usage

monty_dsl_distributions()

Value

A data.frame with columns

• name the name of the distribution; each name begins with a capital letter, and there are duplicate names where different parameterisations are supported.

• args the arguments of all parameters, except the random variable itself which is given as the first argument to density functions.

We may expand the output here in the future to include information on if distributions have support in C++, but we might end up supporting everything this way soon.

Examples

monty_dsl_distributions()

Description

Explain error codes produced by monty. This is a work in progress, and we would like feedback on what is useful as we improve it. The idea is that if you see an error you can link through to get more information on what it means and how to resolve it. The current implementation of this will send you to the rendered vignettes, but in the future we will arrange for offline rendering too.

Usage

monty_dsl_error_explain(code, how = "pretty")

Arguments

Value

Nothing, this is called for its side effect only

Examples

monty_dsl_error_explain("E201")

Description

Parse an expression as if it were a call to one of monty's distribution functions (e.g., Normal, Poisson). This will fill in any defaults, disambiguate where multiple parameterisations of the distribution are available, and provide 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.

Usage

monty_dsl_parse_distribution(expr, name = NULL)

Arguments

Value

A list; the contents of this are subject to change. However you can (to a degree) rely on the following elements:

• name: The name of the distribution (e.g., Normal). This will be the same as the name of the function called in expr

• variant: The name of the distribution variant, if more than one is supported.

• args: The arguments that you provided, in position-matched order

• cpp: The names of the C++ entrypoint to use. This is a list with elements density and sample for the log-density and sampling functions, and NULL where these do not yet exist.

Currently we also include:

• density: A function to compute the log-density. This will likely change once we support creation of differentiable models because we will want to do something with the arguments provided!

• sample: A function to sample from the distribution, given (as a first argument) a rng object (see monty_rng)

Examples

# A successful match
monty_dsl_parse_distribution(quote(Normal(0, 1)))

# An unsuccessful match
monty_dsl_parse_distribution(quote(Normal()))

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.

Usage

monty_example(name, ...)

Arguments

Value

A monty_model object

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.

Examples

monty_example("banana")
monty_example("gaussian", diag(2))

Description

Create a basic monty model. This takes a user-supplied object that minimally can compute a probability density (via a density function) and information about parameters; with this we can sample from the model using MCMC using monty_sample. We don't imagine that many users will call this function directly, but that this will be glue used by packages.

Usage

monty_model(model, properties = NULL)

Arguments

Details

The model argument can be a list or environment (something indexable by $) and have elements:

• density: A function that will compute some probability density. It must take an argument representing a parameter vector (a numeric vector) and return a single value. This is the posterior probability density in Bayesian inference, but it could be anything really. Models can return -Inf if things are impossible, and we'll try and cope gracefully with that wherever possible. If the property allow_multiple_parameters is TRUE, then this function must be able to handle the argument parameter being a matrix, and return a vector of densities.

• parameters: A character vector of parameter names. This vector is the source of truth for the length of the parameter vector.

• domain: Information on the parameter domain. This is a two column matrix with length(parameters) rows representing each parameter. The parameter minimum and maximum bounds are given as the first and second column. Infinite values (-Inf or Inf) should be used where the parameter has infinite domain up or down. Currently used to translate from a bounded to unbounded space for HMC, but we might also use this for reflecting proposals in MCMC too, as well as a fast way of avoiding calculating densities where proposals fall out of bounds. If not present we assume that the model is valid everywhere (i.e., that all parameters are valid from -Inf to Inf. If unnamed, you must provide a domain for all parameters. If named, then you can provide a subset, with parameters that are not included assumed to have a domain of ⁠(-Inf, Inf)⁠.

• direct_sample: A function to sample directly from the parameter space, given a monty_rng object to sample from. In the case where a model returns a posterior (e.g., in Bayesian inference), this is assumed to be sampling from the prior. We'll use this for generating initial conditions for MCMC where those are not given, and possibly other uses. If not given then when using monty_sample() the user will have to provide a vector of initial states.

• gradient: A function to compute the gradient of density with respect to the parameter vector; takes a parameter vector and returns a vector the same length. For efficiency, the model may want to be stateful so that gradients can be efficiently calculated after a density calculation, or density after gradient, where these are called with the same parameters. This function is optional (and may not be well defined or possible to define).

• set_rng_state: A function to set the state (this is in contrast to the rng that is passed through to direct_sample as that is the sampler's rng stream, but we assume models will look after their own stream, and that they may need many streams). Models that provide this method are assumed to be stochastic; however, you can use the is_stochastic property (via monty_model_properties()) to override this (e.g., to run a stochastic model with its deterministic expectation). This function takes a raw vector of random number state from monty_rng and uses it to set the random number state for your model; this is derived from the random number stream for a particular chain, jumped ahead.

• get_rng_state: A function to get the RNG state; must be provided if set_rng_state is present. Must return the random number state, which is a raw vector (potentially quite long).

• parameter_groups: Optionally, an integer vector indicating parameter group membership. The format here may change (especially if we explore more complex nestings) but at present parameters with group 0 affect everything (so are accepted or rejected as a whole), while parameters in groups 1 to n are independent (for example, changing the parameters in group 2 does not affect the density of parameters proposed in group 3).

Value

An object of class monty_model. This will have elements:

• model: The model as provided

• parameters: The parameter name vector

• parameter_groups: The parameter groups

• domain: The parameter domain matrix, named with your parameters

• direct_sample: The direct_sample function, if provided by the model

• gradient: The gradient function, if provided by the model

• properties: A list of properties of the model; see monty_model_properties(). Currently this contains:

  • has_gradient: the model can compute its gradient

  • has_direct_sample: the model can sample from parameters space

  • is_stochastic: the model will behave stochastically

  • has_parameter_groups: The model has separable parameter groups

See Also

monty_model_function, which provides a simple interface for creating models from functions.

Description

Combine two models by multiplication. We'll need a better name here. In Bayesian inference we will want to create a model that represents the multiplication of a likelihood and a prior (in log space) and it will be convenient to think about these models separately. Multiplying probabilities (or adding on a log scale) is common enough that there may be other situations where we want to do this.

Usage

monty_model_combine(a, b, properties = NULL, name_a = "a", name_b = "b")

Arguments

Details

Here we describe the impact of combining a pair of models

• density: this is the sum of the log densities from each model

• parameters: the union of parameters from each model is taken

• domain: The most restrictive domain is taken for each parameter. Parameters that do not appear in one model are assumed to have infinite domain there.

• gradient: if both models define a gradient, this is the sum of the gradients. If either does not define a gradient, the resulting model will not have gradient support. Set has_gradient = TRUE within 'properties if you want to enforce that the combination is differentiable. If the models disagree in their parameters, parameters that are missing from a model are assumed (reasonably) to have a zero gradient.

• direct_sample: this one is hard to do the right thing for. If neither model can be directly sampled from that's fine, we don't directly sample. If only one model can be sampled from and if it can sample from the union of all parameters then we take that function (this is the case for a prior model when combined with a likelihood). Other cases will be errors, which can be avoided by setting has_direct_gradient = FALSE in properties.

• is_stochastic: a model is stochastic if either component is stochastic.

The properties of the model will be combined as above, reflecting the properties of the joint model.

The model field will be an ordered, unnamed, list containing the two elements corresponding to the first and second model (not the monty_model, but the underlying model, perhaps?). This is the only part that makes a distinction between the two models here; for all components above they are equivalent.

Value

A monty_model object

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))

Description

Compute log density for a model. This is a wrapper around the ⁠$density⁠ property within a monty_model object.

Usage

monty_model_density(model, parameters)

Arguments

Value

A log-density value, or vector of log-density values

See Also

monty_model_gradient for computing gradients and monty_model_direct_sample for sampling from a model.

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))

Description

Directly sample from a model. Not all models support this, and an error will be thrown if it is not possible.

Usage

monty_model_direct_sample(model, rng, named = FALSE)

Arguments

Value

A vector or matrix of sampled parameters

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)

Description

Create a monty_model from a function that computes density. This allows use of any R function as a simple monty model. If you need advanced model features, then this interface may not suit you and you may prefer to use monty_model directly.

Usage

monty_model_function(
  density,
  packer = NULL,
  fixed = NULL,
  allow_multiple_parameters = FALSE
)

Arguments

Details

This interface will expand in future versions of monty to support gradients, stochastic models, parameter groups and simultaneous calculation of density.

Value

A monty_model object that computes log density with the provided density function, given a numeric vector argument representing all parameters.

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))

Description

Compute the gradient of log density (which is returned by monty_model_density) with respect to parameters. Not all models support this, and an error will be thrown if it is not possible.

Usage

monty_model_gradient(model, parameters, named = FALSE)

Arguments

Value

A vector or matrix of gradients

See Also

monty_model_density for log density, and monty_model_direct_sample to sample from a model

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))

Description

Describe properties of a model. Use of this function is optional, but you can pass the return value of this as the properties argument of monty_model to enforce that your model does actually have these properties.

Usage

monty_model_properties(
  has_gradient = NULL,
  has_direct_sample = NULL,
  is_stochastic = NULL,
  has_parameter_groups = NULL,
  has_observer = NULL,
  allow_multiple_parameters = NULL
)

Arguments

Value

A list of class monty_model_properties which should not be modified.

Examples

# Default properties:
monty_model_properties()

# Set some properties:
monty_model_properties(has_gradient = TRUE, is_stochastic = FALSE)

Description

Split a model that has been combined by monty_model_combine() into its constituent parts.

Usage

monty_model_split(model, prior_first = FALSE)

Arguments

Details

We assume that a split model can be broken into a "prior" and a "likelihood" if exactly one model:

• can be directly sampled from

• is not stochastic

• consumes all parameters

Typically, it will be the first criterion that will separate a model into prior and likelihood (if you could sample from your likelihood, then you would not use a sampler, which is where we are typically going to perform this action).

If prior_first is FALSE we just return the parts.

Value

An unnamed list of length 2, being the component models. If one model might be the prior it will be listed first.

Description

Create an observer to extract additional details from your model during the sampling process.

Usage

monty_observer(observe, finalise = NULL, combine = NULL, append = NULL)

Arguments

Details

Sometimes you want to extract additional information from your model as your chain runs. The case we see this most is when running MCMC with a particle filter (pmcmc); in this case while the likelihood calculation is running we are computing lots of interesting quantities such as the final state of the system (required for onward simulation) and filtered trajectories through time. Because these are stochastic we can't even just rerun the model with our sampled parameter sets, because the final states that are recovered depend also on the random number generators (practically we would not want to, as it is quite expensive to compute these quantities).

The observer mechanism allows you to carry out arbitrary additional calculations with your model at the end of the step.

Value

An object with class monty_observer which can be passed in to monty_sample.

Description

Build a packer, which can be used to translate between an unstructured vector of numbers (the vector being updated by an MCMC for example) and a structured list of named values, which is easier to program against. This is useful for the bridge between model parameters and a model's implementation, but it is also useful for the state vector in a state-space model. We refer to the process of taking a named list of scalars, vectors and arrays and converting into a single vector "packing" and the inverse "unpacking".

Usage

monty_packer(scalar = NULL, array = NULL, fixed = NULL, process = NULL)

Arguments

Details

There are several places where it is most convenient to work in an unstructured vector:

• An MCMC is typically discussed as a the updating of some vector x to another ⁠x'⁠

• An optimisation algorithm will try and find a set of values for a vector x that minimises (or maximises) some function f(x)

• An ode solver works with a vector x(t) (x at time t) and considers x(t + h) by computing the vector of derivatives dx(t)/dt

In all these cases, the algorithm that needs the vector of numbers knows nothing about what they represent. Commonly, these will be a packed vector of parameters. So our vector x might actually represent the parameters a, b and c in a vector as ⁠[a, b, c]⁠ - this is a very common pattern, and you have probably implemented this yourself.

In more complex settings, we might want our vector x to collect more structured quantities. Suppose that you are fitting a model with an age-structured or sex-structured parameter. Rather than having a series of scalars packed into your vector x you might have a series of values destined to be treated as a vector:

| 1  2  3  4  5  6  7  |
| a  b  c  d1 d2 d3 d4 |

So here we might have a vector of length 7, where the first three elements will represent be the scalar values a, b and c but the next four will be a vector d.

Unpacked, this might be written as:

list(a = 1, b = 2, c = 3, d = 4:7)

The machinery here is designed to make these transformations simple and standardised within monty, and should be flexible enough for many situations. We will also use these from within dust2 and odin2 for transformations in and out of vectors of ODE state.

Value

An object of class monty_packer, which has elements:

• names: a function that returns a character vector of computed names; in the parameter packer context these are the names that your statistical model will use.

• unpack: a function that can unpack an unstructured vector (say, from your statistical model parameters) into a structured list (say, for your generative model)

• pack: a function that can pack your structured list of data back into a numeric vector, for example suitable for a statistical model. This ignores values created by a preprocess function and present in fixed.

• index: a function which produces a named list where each element has the name of a value in scalar or array and each value has the indices within an unstructured vector where these values can be found, in the shape of the data that would be unpacked. This is of limited most use to most people.

• subset: an experimental interface which can be used to subset a packer to a packer for a subset of contents. Documentation will be provided once the interface settles, but this is for advanced use only!

When to use process

The process function is a get-out-of-jail function designed to let you do arbitrary transformations when unpacking a vector. In general, this should not be the first choice to use because it is less easy to reason about by other tooling (for example, as we develop automatic differentiation support for use with the HMC algorithm, a process function will be problematic because we will need to make sure we can differentiate this process). However, there are cases where it will be only way to achieve some results.

Imagine that you are packing a 2x2 covariance matrix into your vector in order to use within an MCMC or optimisation algorithm. 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 array, so that our unpacked list holds b_raw = c(b11, b12, b22). We would then write process as something like:

process <- function(x) {
  list(b = matrix(x$b_raw[c(1, 2, 2, 3)], 2, 2))
}

which creates the symmetric 2x2 matrix b from b_raw.

Unpacking matrices

If you do not use fixed or process when defining your packer, then you can use ⁠$unpack()⁠ with a matrix or higher-dimensional output. There are two ways that you might like to unpack this sort of output. Assume you have a matrix m with 3 rows and 2 columns; this means that we have two sets of parameters or state (one per column) and 3 states within each; this is the format that MCMC parameters will be in for example.

The first would to be return a list where the ith element is the result of unpacking the ith parameter/state vector. You can do this by running

apply(m, 2, p$unpack)

The second would be to return a named list with three elements where the ith element is the unpacked version of the ith state. In this case you can pass the matrix directly in to the unpacker:

p$unpack(m)

When you do this, the elements of m will acquire an additional dimension; scalars become vectors (one per set), vectors become matrices (one column per set) and so on.

This approach generalises to higher dimensional input, though we suspect you'll spend a bit of time head-scratching if you use it.

Packing lists into vectors and matrices

The unpacking operation is very common - an MCMC proceeds, produces an unstructured vector, and you unpack it into a list in order to be able to easily work with it. The reverse is much less common, where we take a list and convert it into a vector (or matrix, or multidimensional array). Use of this direction ("packing") may be more common where using packers to work with the output of state-space models (e.g. in odin2 or dust2, which use this machinery).

The input to pack() will be the shape that unpack() returned; a named list of numerical vectors, matrices and arrays. The names must correspond to the names in your packer (i.e., scalar and the names of array). Each element has dimensions

<...object, ...residual>

where ...object is the dimensions of the data itself and ...residual is the dimensions of the hypothetical input to pack.

There is an unfortunate ambiguity in R's lack of true scalar types that we cannot avoid. It is hard to tell the difference packing a vector vs packing an array where all dimensions are 1. See the examples, and please let us know if the behaviour needs changing.

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))

# R lacks scalars, which means that some packers will unpack
# different inputs to the same outputs:
p <- monty_packer(c("a", "b"))
p$unpack(1:2)
p$unpack(cbind(1:2))

# This means that we can't reliably pack these inputs in a way
# that guarantees round-tripping is possible.  We have chosen to
# prioritise the case where a *single vector* is round-trippable:
p$pack(list(a = 1, b = 2))

# This ambiguity goes away if unpacking matices with more than one
# column:
p$unpack(matrix(1:6, 2, 3))

Description

Build a grouped version of monty_packer() with the same basic idea; convert between a vector representation of some group of numbers to a named list of structured data, but with an extra twist: this time the unstructured vector of numbers contains values that correspond to multiple groups and some are shared across groups while others vary between groups. This function does a lot of bookkeeping in a relatively short amount of code, so you should be familiar with the ideas in monty_packer() before continuing.

Usage

monty_packer_grouped(
  groups,
  scalar = NULL,
  array = NULL,
  fixed = NULL,
  process = NULL,
  shared = NULL
)

Arguments

Details

Recall from monty_packer() that our original problem was to take an unstructured vector like

| 1  2  3  4  5  6  7  |
| a  b  c  d1 d2 d3 d4 |

and unpack it into a structured list like

list(a = 1, b = 2, c = 3, d = 4:7)

Our aim here is to do the same but to allow some of these values (say b and c) to be shared (constant) over groups while the others (a and d) to vary by group. So for groups x and y we might try and create something like

list(
  list(x = list(a = 3, b = 1, c = 2, d = 4:7),
       y = list(a = 8, b = 1, c = 2, d = 9:12))

from a vector

| 1  2  3  4  5  6  7  8  9  10 11 12 |
| b  c  a  d1 d2 d3 d4 a  d1 d2 d3 d4 |
| xy xy x  x  x  x  x  y  y  y  y  y  |

Value

An object of class monty_packer_grouped, which has the same elements as monty_packer, though with slightly different effects.

• names: a function that returns a character vector of computed names; in the parameter packer context these are the names that your statistical model will use.

• groups: A function that returns your group names (the groups argument as supplied)

• unpack: A function for converting from an unstructured vector into a nested list. Each element of this list is conceptually the same as the result of unpack() from monty_packer().

• pack: The inverse to unpack() but less commonly performed. Convert a nested list into an unstructured vector. Quite a lot of validation is required to make sure that the input has not been tampered with, and errors thrown while doing this validation may not be very interpretable.

• index: The nested version of the index() function in monty_packer(). The outer list is over groups, and the inner list contains the position within the original unstructured vector where this value can be found. It is not clear to us if this is a useful list.

• subset: A function that might eventually allow subsetting a grouped packer. Currently it just errors.

Examples

p <- monty_packer_grouped(c("x", "y"), c("a", "b", "c", "d", "e"),
                          shared = c("b", "c"))
p$names()
p$unpack(1:8)

Description

Create an object that can be used to generate random numbers with the same RNG as monty uses internally. This is primarily meant for debugging and testing the underlying C++ rather than a source of random numbers from R.

Value

A monty_rng object, which can be used to drawn random numbers from monty's distributions.

Warning

This interface is subject to change in the near future, we do not recommend its use in user code.

Running multiple streams, perhaps in parallel

The underlying random number generators are designed to work in parallel, and with random access to parameters (see vignette("rng") for more details). However, this is usually done within the context of running a model where each particle sees its own stream of numbers. We provide some support for running random number generators in parallel, but any speed gains from parallelisation are likely to be somewhat eroded by the overhead of copying around a large number of random numbers.

All the random distribution functions support an argument n_threads which controls the number of threads used. This argument will silently have no effect if your installation does not support OpenMP.

Parallelisation will be performed at the level of the stream, where we draw n numbers from each stream for a total of n * n_streams random numbers using n_threads threads to do this. Setting n_threads to be higher than n_streams will therefore have no effect. If running on somebody else's system (e.g., an HPC, CRAN) you must respect the various environment variables that control the maximum allowable number of threads.

With the exception of random_real, each random number distribution accepts parameters; the interpretations of these will depend on n, n_streams and their rank.

• If a scalar then we will use the same parameter value for every draw from every stream

• If a vector with length n then we will draw n random numbers per stream, and every stream will use the same parameter value for every stream for each draw (but a different, shared, parameter value for subsequent draws).

• If a matrix is provided with one row and n_streams columns then we use different parameters for each stream, but the same parameter for each draw.

• If a matrix is provided with n rows and n_streams columns then we use a parameter value ⁠[i, j]⁠ for the ith draw on the jth stream.

The rules are slightly different for the prob argument to multinomial as for that prob is a vector of values. As such we shift all dimensions by one:

• If a vector we use same prob every draw from every stream and there are length(prob) possible outcomes.

• If a matrix with n columns then vary over each draw (the ith draw using vector prob[, i] but shared across all streams. There are nrow(prob) possible outcomes.

• If a 3d array is provided with 1 column and n_streams "layers" (the third dimension) then we use then we use different parameters for each stream, but the same parameter for each draw.

• If a 3d array is provided with n columns and n_streams "layers" then we vary over both draws and streams so that with use vector prob[, i, j] for the ith draw on the jth stream.

The output will not differ based on the number of threads used, only on the number of streams.

Public fields

Methods

Public methods

• monty_rng$new()

• monty_rng$size()

• monty_rng$jump()

• monty_rng$long_jump()

• monty_rng$random_real()

• monty_rng$random_normal()

• monty_rng$uniform()

• monty_rng$normal()

• monty_rng$binomial()

• monty_rng$beta_binomial_ab()

• monty_rng$beta_binomial_prob()

• monty_rng$negative_binomial_prob()

• monty_rng$negative_binomial_mu()

• monty_rng$hypergeometric()

• monty_rng$gamma_scale()

• monty_rng$gamma_rate()

• monty_rng$poisson()

• monty_rng$exponential_rate()

• monty_rng$exponential_mean()

• monty_rng$cauchy()

• monty_rng$multinomial()

• monty_rng$beta()

• monty_rng$truncated_normal()

• monty_rng$state()

• monty_rng$set_state()

Method new()

Create a monty_rng object

Usage

monty_rng$new(seed = NULL, n_streams = 1L, deterministic = FALSE)

Arguments

Method size()

Number of streams available

Usage

monty_rng$size()

Method jump()

The jump function updates the random number state for each stream by advancing it to a state equivalent to 2^128 numbers drawn from each stream.

Usage

monty_rng$jump()

Method long_jump()

Longer than ⁠$jump⁠, the ⁠$long_jump⁠ method is equivalent to 2^192 numbers drawn from each stream.

Usage

monty_rng$long_jump()

Method random_real()

Generate n numbers from a standard uniform distribution

Usage

monty_rng$random_real(n, n_threads = 1L)

Arguments

Method random_normal()

Generate n numbers from a standard normal distribution

Usage

monty_rng$random_normal(n, n_threads = 1L, algorithm = "box_muller")

Arguments

Method uniform()

Generate n numbers from a uniform distribution

Usage

monty_rng$uniform(n, min, max, n_threads = 1L)

Arguments

Method normal()

Generate n numbers from a normal distribution

Usage

monty_rng$normal(n, mean, sd, n_threads = 1L, algorithm = "box_muller")

Arguments

Method binomial()

Generate n numbers from a binomial distribution

Usage

monty_rng$binomial(n, size, prob, n_threads = 1L)

Arguments

Method beta_binomial_ab()

Generate n numbers from a beta-binomial distribution

Usage

monty_rng$beta_binomial_ab(n, size, a, b, n_threads = 1L)

Arguments

Method beta_binomial_prob()

Generate n numbers from a beta-binomial distribution

Usage

monty_rng$beta_binomial_prob(n, size, prob, rho, n_threads = 1L)

Arguments

Method negative_binomial_prob()

Generate n numbers from a negative binomial distribution

Usage

monty_rng$negative_binomial_prob(n, size, prob, n_threads = 1L)

Arguments

Method negative_binomial_mu()

Generate n numbers from a negative binomial distribution

Usage

monty_rng$negative_binomial_mu(n, size, mu, n_threads = 1L)

Arguments

Method hypergeometric()

Generate n numbers from a hypergeometric distribution

Usage

monty_rng$hypergeometric(n, n1, n2, k, n_threads = 1L)

Arguments

Method gamma_scale()

Generate n numbers from a gamma distribution

Usage

monty_rng$gamma_scale(n, shape, scale, n_threads = 1L)

Arguments

Method gamma_rate()

Generate n numbers from a gamma distribution

Usage

monty_rng$gamma_rate(n, shape, rate, n_threads = 1L)

Arguments

Method poisson()

Generate n numbers from a Poisson distribution

Usage

monty_rng$poisson(n, lambda, n_threads = 1L)

Arguments

Method exponential_rate()

Generate n numbers from a exponential distribution

Usage

monty_rng$exponential_rate(n, rate, n_threads = 1L)

Arguments

Method exponential_mean()

Generate n numbers from a exponential distribution

Usage

monty_rng$exponential_mean(n, mean, n_threads = 1L)

Arguments

Method cauchy()

Generate n draws from a Cauchy distribution.

Usage

monty_rng$cauchy(n, location, scale, n_threads = 1L)

Arguments

Method multinomial()

Generate n draws from a multinomial distribution. In contrast with most of the distributions here, each draw is a vector with the same length as prob.

Usage

monty_rng$multinomial(n, size, prob, n_threads = 1L)

Arguments

Method beta()

Generate n numbers from a beta distribution

Usage

monty_rng$beta(n, a, b, n_threads = 1L)

Arguments

Method truncated_normal()

Generate n numbers from a truncated normal distribution

Usage

monty_rng$truncated_normal(n, mean, sd, min, max, n_threads = 1L)

Arguments

Method state()

Returns the state of the random number stream. This returns a raw vector of length 32 * n_streams.

Usage

monty_rng$state()

Method set_state()

Sets the state of the random number stream.

Usage

monty_rng$set_state(state)

Arguments

Examples

rng <- monty::monty_rng$new(42)

# Shorthand for Uniform(0, 1)
rng$random_real(5)

# Shorthand for Normal(0, 1)
rng$random_normal(5)

# Uniform random numbers between min and max
rng$uniform(5, -2, 6)

# Normally distributed random numbers with mean and sd
rng$normal(5, 4, 2)

# Binomially distributed random numbers with size and prob
rng$binomial(5, 10, 0.3)

# Negative binomially distributed random numbers with size and prob
rng$negative_binomial_prob(5, 10, 0.3)

# Negative binomially distributed random numbers with size and mean mu
rng$negative_binomial_mu(5, 10, 25)

# Hypergeometric distributed random numbers with parameters n1, n2 and k
rng$hypergeometric(5, 6, 10, 4)

# Gamma distributed random numbers with parameters shape and scale
rng$gamma_scale(5, 0.5, 2)

# Gamma distributed random numbers with parameters shape and rate
rng$gamma_rate(5, 0.5, 2)

# Poisson distributed random numbers with mean lambda
rng$poisson(5, 2)

# Exponentially distributed random numbers with rate
rng$exponential_rate(5, 2)

# Exponentially distributed random numbers with mean
rng$exponential_mean(5, 0.5)

# Multinomial distributed random numbers with size and vector of
# probabiltiies prob
rng$multinomial(5, 10, c(0.1, 0.3, 0.5, 0.1))

Description

Create a set of initial random number seeds suitable for using within a distributed context (over multiple processes or nodes) at a level higher than a single group of synchronised threads.

Usage

monty_rng_distributed_state(
  seed = NULL,
  n_streams = 1L,
  n_nodes = 1L,
  algorithm = "xoshiro256plus"
)

monty_rng_distributed_pointer(
  seed = NULL,
  n_streams = 1L,
  n_nodes = 1L,
  algorithm = "xoshiro256plus"
)

Arguments

Details

See vignette("rng_distributed") for a proper introduction to these functions.

Value

A list of either raw vectors (for monty_rng_distributed_state) or of monty_rng_pointer objects (for monty_rng_distributed_pointer)

Examples

monty::monty_rng_distributed_state(n_nodes = 2)
monty::monty_rng_distributed_pointer(n_nodes = 2)

Description

This function exists to support use from other packages that wish to use monty's random number support, and creates an opaque pointer to a set of random number streams.

Public fields

Methods

Public methods

• monty_rng_pointer$new()

• monty_rng_pointer$sync()

• monty_rng_pointer$state()

• monty_rng_pointer$is_current()

Method new()

Create a new monty_rng_pointer object

Usage

monty_rng_pointer$new(
  seed = NULL,
  n_streams = 1L,
  long_jump = 0L,
  algorithm = "xoshiro256plus"
)

Arguments

Method sync()

Synchronise the R copy of the random number state. Typically this is only needed before serialisation if you have ever used the object.

Usage

monty_rng_pointer$sync()

Method state()

Return a raw vector of state. This can be used to create other generators with the same state.

Usage

monty_rng_pointer$state()

Method is_current()

Return a logical, indicating if the random number state that would be returned by state() is "current" (i.e., the same as the copy held in the pointer) or not. This is TRUE on creation or immediately after calling ⁠$sync()⁠ or ⁠$state()⁠ and FALSE after any use of the pointer.

Usage

monty_rng_pointer$is_current()

Examples

monty::monty_rng_pointer$new()

Description

Run MCMC chains in parallel (at the same time). This runner uses the callr package to distribute your chains over a number of worker processes on the same machine. If you have used mcstate, this is the same as "worker" processes. Unless your chains take a few seconds to run, this will be slower than running with the default serial runner (monty_runner_serial), however for long running chains, the speed-up will typically scale with workers added, so long as your chains divide neatly over workers.

Usage

monty_runner_callr(n_workers, progress = NULL)

Arguments

Value

A runner of class monty_runner that can be passed to monty_sample()

Description

Run MCMC chains in parallel (at the same time). This runner uses the parallel package to distribute your chains over a number of worker processes on the same machine. Compared with monty_runner_callr (Which is similar to the "worker" support in mcstate version 1), this is very simple. In particular we do not report back any information about progress while a chain is running on a worker or even across chains. There's also no support to warn you if your number of chains do not neatly divide through by the number of workers. Mostly this exists as a proof of concept for us to think about the different interfaces. Unless your chains are quite slow, the parallel runner will be slower than the serial runner (monty_runner_serial) due to the overhead cost of starting the cluster.

Usage

monty_runner_parallel(n_workers)

Arguments

Value

A runner of class monty_runner that can be passed to monty_sample()

Description

Run MCMC chains in series (one after another). This is the simplest chain runner, and the default used by monty_sample(). It has nothing that can be configured (yet).

Usage

monty_runner_serial(progress = NULL)

Arguments

Value

A runner of class monty_runner that can be passed to monty_sample()

Description

Run chains simultaneously. This differs from monty_runner_parallel, which runs chains individually in parallel by working with models that can evaluate multiple densities at the same time. There are situations where this might be faster than running in parallel, but primarily this exists so that we can see that samplers can work with multiple samples at once.

Usage

monty_runner_simultaneous(progress = NULL)

Arguments

Value

A runner of class monty_runner that can be passed to monty_sample()

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)

Description

Sample from a model. Uses a Monte Carlo method (or possibly something else in future) to generate samples from your distribution. This is going to change a lot in future, as we add 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.

Usage

monty_sample(
  model,
  sampler,
  n_steps,
  initial = NULL,
  n_chains = 1L,
  runner = NULL,
  restartable = FALSE,
  burnin = NULL,
  thinning_factor = NULL
)

Arguments

Value

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 ith parameter from the jth sample from the kth chain. The rows will be named with the names of the parameters, from your model.

• density: A matrix of model log densities, with n_steps rows and n_chains columns.

• initial: A record of the initial conditions, a matrix with as many rows as you have parameters and n_chains columns (this is the same format as the matrix form of the initial input parameter)

• details: Additional details reported by the sampler; this will be a list of length n_chains (or NULL) and the details depend on the sampler. This one is subject to change.

• observations: Additional details reported by the model. This one is also subject to change.

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:

res <- posterior::as_draws_df(samples)
posterior::summarise_draws(res)

# At this point you could also use the 'bayesplot' package to plot
# diagnostics.

Description

Continue (restart) chains started by monty_sample. Requires that the original chains were run with restartable = TRUE. Running chains this way will result in the final state being exactly the same as running for the total (original + continued) number of steps in a single push.

Usage

monty_sample_continue(
  samples,
  n_steps,
  restartable = FALSE,
  runner = NULL,
  append = TRUE
)

Arguments

Value

A list of parameters and densities

Description

Clean up after manual sampling. This is essentially a safe version of deleting the directory (e.g, unlink(path, recursive = TRUE)) which checks that the directory really was used for sampling and that it does not contain anything else unexpected.

Usage

monty_sample_manual_cleanup(path)

Arguments

Value

Nothing, called for side effects only.

Examples

model <- monty_example("banana")
sampler <- monty_sampler_random_walk(vcv = diag(2) * 0.05)
path <- tempfile()
monty_sample_manual_prepare(model, sampler, 100, path)
monty_sample_manual_info(path)

# Run the (single) chain
monty_sample_manual_run(1, path)
monty_sample_manual_info(path)

# Collect the results
monty_sample_manual_collect(path)

# Clean up samples
monty_sample_manual_cleanup(path)

Description

Collect samples from chains that have been run with monty_sample_manual_prepare and monty_sample_manual_run. If any chain has not completed, we will error.

Usage

monty_sample_manual_collect(
  path,
  samples = NULL,
  restartable = FALSE,
  append = TRUE
)

Arguments

Value

A monty_samples object.

Examples

model <- monty_example("banana")
sampler <- monty_sampler_random_walk(vcv = diag(2) * 0.05)
path <- tempfile()
monty_sample_manual_prepare(model, sampler, 100, path)
monty_sample_manual_info(path)

# Run the (single) chain
monty_sample_manual_run(1, path)
monty_sample_manual_info(path)

# Collect the results
monty_sample_manual_collect(path)

# Clean up samples
monty_sample_manual_cleanup(path)

Description

Get information about the status of manually scheduled samples.

Usage

monty_sample_manual_info(path)

Arguments

Value

Invisibly, a logical vector, TRUE for completed chains and FALSE for incomplete chains.

Examples

model <- monty_example("banana")
sampler <- monty_sampler_random_walk(vcv = diag(2) * 0.05)
path <- tempfile()
monty_sample_manual_prepare(model, sampler, 100, path)
monty_sample_manual_info(path)

# Run the (single) chain
monty_sample_manual_run(1, path)
monty_sample_manual_info(path)

# Collect the results
monty_sample_manual_collect(path)

# Clean up samples
monty_sample_manual_cleanup(path)

Description

Run an MCMC, but schedule execution of the chains yourself. Use this if you want to distribute chains over (say) the nodes of an HPC system. The arguments are the same as for monty_sample, except that the runner argument is missing as you will be looking after that yourself. After using this function, you will generally be wanting to run monty_sample_manual_run and monty_sample_manual_collect.

Usage

monty_sample_manual_prepare(
  model,
  sampler,
  n_steps,
  path,
  initial = NULL,
  n_chains = 1L,
  burnin = NULL,
  thinning_factor = NULL
)

Arguments

Details

In contrast to monty_sample there is no runner argument here, because by using this function directly you are taking responsibility for being your own runner.

As with the ways of running a set of chains in monty, it is expected that using ⁠monty_sample_manual_*⁠ will result in the same samples being generated as if you had used monty_sample with a runner of your choice.

Value

Invisibly, the path used to store files (the same as the value of the path argument)

See Also

monty_sample_manual_run to run the chains and monty_sample_manual_collect / monty_sample_manual_cleanup to collect results and clean up. monty_sample_manual_info can print human-readable information about the state of a manual run.

Examples

model <- monty_example("banana")
sampler <- monty_sampler_random_walk(vcv = diag(2) * 0.05)
path <- tempfile()
monty_sample_manual_prepare(model, sampler, 100, path)
monty_sample_manual_info(path)

# Run the (single) chain
monty_sample_manual_run(1, path)
monty_sample_manual_info(path)

# Collect the results
monty_sample_manual_collect(path)

# Clean up samples
monty_sample_manual_cleanup(path)

Description

Prepare to continue sampling from a model, with manual chain orchestration. This function is to monty_sample_continue what monty_sample_manual_prepare is to monty_sample. The original set of samples do not need to have been run manually.

Usage

monty_sample_manual_prepare_continue(
  samples,
  n_steps,
  path,
  save_samples = "hash"
)

Arguments

Value

Invisibly, the path used to store files (the same as the value of the path argument)

Description

Run a chain that was prepared using monty_sample_manual_prepare.

Usage

monty_sample_manual_run(chain_id, path, progress = NULL)

Arguments

Warning:

There is no lock mechanism; you can start a single chain many times. Don't do this.

Examples

model <- monty_example("banana")
sampler <- monty_sampler_random_walk(vcv = diag(2) * 0.05)
path <- tempfile()
monty_sample_manual_prepare(model, sampler, 100, path)
monty_sample_manual_info(path)

# Run the (single) chain
monty_sample_manual_run(1, path)
monty_sample_manual_info(path)

# Collect the results
monty_sample_manual_collect(path)

# Clean up samples
monty_sample_manual_cleanup(path)

Description

Create an adaptive Metropolis-Hastings sampler, which will tune its variance covariance matrix (vs the simple random walk sampler monty_sampler_random_walk).

Usage

monty_sampler_adaptive(
  initial_vcv,
  initial_vcv_weight = 1000,
  initial_scaling = 1,
  initial_scaling_weight = NULL,
  min_scaling = 0,
  scaling_increment = NULL,
  log_scaling_update = TRUE,
  acceptance_target = 0.234,
  forget_rate = 0.2,
  forget_end = Inf,
  adapt_end = Inf,
  pre_diminish = 0,
  boundaries = "reflect"
)

Arguments

Details

Efficient exploration of the parameter space during an MCMC might be difficult when the target distribution is of high dimensionality, especially if the target probability distribution present a high degree of correlation. Adaptive schemes are used to "learn" on the fly the correlation structure by updating the proposal distribution by recalculating the empirical variance-covariance matrix and rescale it at each adaptive step of the MCMC.

Our implementation of an adaptive MCMC algorithm is based on an adaptation of the "accelerated shaping" algorithm in Spencer (2021). The algorithm is based on a random-walk Metropolis-Hastings algorithm where the proposal is a multi-variate Normal distribution centred on the current point.

Spencer SEF (2021) Accelerating adaptation in the adaptive Metropolis–Hastings random walk algorithm. Australian & New Zealand Journal of Statistics 63:468-484.

Value

A monty_sampler object, which can be used with monty_sample

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


coda::effectiveSize(coda::as.mcmc.list(res_rw))
coda::effectiveSize(coda::as.mcmc.list(res_adapt))

Description

Create a Hamiltonian Monte Carlo sampler, implemented using the leapfrog algorithm.

Usage

monty_sampler_hmc(
  epsilon = 0.015,
  n_integration_steps = 10,
  vcv = NULL,
  debug = FALSE
)

Arguments

Value

A monty_sampler object, which can be used with monty_sample

Description

Create a nested adaptive Metropolis-Hastings sampler, which extends the adaptive sampler monty_sampler_adaptive, tuning the variance covariance matrices for proposal for the separable sections of a nested model (vs the simple nested random walk sampler monty_sampler_random_walk). This sampler requires that models support the has_parameter_groups property.

Usage

monty_sampler_nested_adaptive(
  initial_vcv,
  initial_vcv_weight = 1000,
  initial_scaling = 1,
  initial_scaling_weight = NULL,
  min_scaling = 0,
  scaling_increment = NULL,
  log_scaling_update = TRUE,
  acceptance_target = 0.234,
  forget_rate = 0.2,
  forget_end = Inf,
  adapt_end = Inf,
  pre_diminish = 0,
  boundaries = "reflect"
)

Arguments

Details

Much like the simple nested random walk sampler monty_sampler_random_walk, the strategy is to propose all the shared parameters as a deviation from the current point in parameter space as a single move and accept or reject as a block. Then we generate points for all the region-specific parameters, compute the density and then accept or reject these updates independently. This is possible because the change in likelihood in region A is independent from region B.

The adaptive proposal algorithm of the non-nested adaptive sampler monty_sampler_adaptive is extended here to adaptively tune the variance covariance matrix of each of these parameter chunks.

Value

A monty_sampler object, which can be used with monty_sample

Description

Create a nested random walk sampler, which uses a symmetric proposal for separable sections of a model to move around in parameter space. This sampler supports sampling from models where the likelihood is only computable randomly (e.g., for pmcmc), and requires that models support the has_parameter_groups property.

Usage

monty_sampler_nested_random_walk(vcv, boundaries = "reflect")

Arguments

Details

The intended use case for this sampler is for models where the density can be decomposed at least partially into chunks that are independent from each other. Our motivating example for this is a model of COVID-19 transmission where some parameters region-specific (e.g., patterns and rates of contact between individuals), and some parameters are shared across all regions (e.g., intrinsic properties of the disease such as incubation period).

The strategy is to propose all the shared parameters as a deviation from the current point in parameter space as a single move and accept or reject as a block. Then we generate points for all the region-specific parameters, compute the density and then accept or reject these updates independently. This is possible because the change in likelihood in region A is independent from region B.

We expect that this approach will be beneficial in limited situations, but where it is beneficial it is likely to result in fairly large speed-ups:

• You probably need more than three regions; as the number of regions increases the benefit of independently accepting or rejecting densities increases (with 1000 separate regions your chains will mix very slowly for example).

• Your model is fairly computationally heavy so that the density calculation completely dominates the sampling process.

• You do not have access to gradient information for your model; we suspect that HMC will outperform this approach by some margin because it already includes this independence via the gradients.

• You can compute your independent calculations in parallel, which help this method reduce your walk time.

Value

A monty_sampler object, which can be used with monty_sample

Description

Create a simple random walk sampler, which uses a symmetric proposal to move around parameter space. This sampler supports sampling from models where the likelihood is only computable randomly (e.g., for pmcmc).

Usage

monty_sampler_random_walk(
  vcv = NULL,
  boundaries = "reflect",
  rerun_every = Inf,
  rerun_random = TRUE
)

Arguments

Value

A monty_sampler object, which can be used with monty_sample

Description

Thin results of running monty_sample(), reducing autocorrelation between samples and saving space. This function may be useful before running onward simulations, or before saving output to disk.

Usage

monty_samples_thin(samples, thinning_factor = NULL, burnin = NULL)

Arguments

Value

A monty_samples object (as for monty_sample()), typically with fewer samples.

Limitations

Subsetting parameters (⁠$pars⁠) and density (⁠$density⁠) is easy enough, but the main use of this function is subsetting chains that have observations, otherwise you could simply cast to samples_df and use functions from the posterior package.

We can only subset observations where the observer was able to tidy them up into a nice array for us. This will typically be the case (for example when using odin/dust, trajectories will be in a nice array).

More specifically, an array is "nice" if the last two dimensions represent samples and chains; in that case we subset along the samples dimension and leave everything else alone. For each element of ⁠$observations⁠ that cannot be subsetted, we will issue a warning.

We cannot generally subset "details", and will pass that along unmodified.

Description

Trace calls to R's random-number-generating functions, to detect unexpected use of random number generation outside of monty's control.

Usage

with_trace_random(code, max_calls = 5, show_stack = FALSE)

Arguments

Value

The result of evaluating code