| Title: | Monte Carlo Models |
|---|---|
| Description: | Experimental sources for the next generation of mcstate, now called 'monty', which will support much of the old mcstate functionality but new things like better parameter interfaces, Hamiltonian Monte Carlo, and other features. |
| Authors: | Rich FitzJohn [aut, cre], Wes Hinsley [aut], Ed Knock [aut], Marc Baguelin [aut], Imperial College of Science, Technology and Medicine [cph] |
| Maintainer: | Rich FitzJohn <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.22 |
| Built: | 2024-12-16 10:24:08 UTC |
| Source: | https://github.com/mrc-ide/monty |
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.
monty_differentiation()monty_differentiation()
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.
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.
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.
We may need to make this slightly extensible in future, but for now the set of functions that can be differentiated is closed.
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.
d <- monty_differentiation() d$differentiate(quote(sqrt(sin(x))), "x") D(quote(sqrt(sin(x))), "x")d <- monty_differentiation() d$differentiate(quote(sqrt(sin(x))), "x") D(quote(sqrt(sin(x))), "x")
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.
monty_domain_expand(domain, packer)monty_domain_expand(domain, packer)
domain |
A two-column matrix as defined in monty_model,
with row names corresponding to either logical names (e.g., |
packer |
A monty_packer object |
A two dimensional matrix representing your domain, or
NULL if domain was given as NULL.
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)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)
Create a model using the monty DSL; this function will likely change name in future, as will its interface.
monty_dsl(x, type = NULL, gradient = NULL, fixed = NULL)monty_dsl(x, type = NULL, gradient = NULL, fixed = NULL)
x |
The model as an expression. This may be given as an
expression, as a string, or as a path to a filename. Typically,
we'll do a reasonable job of working out what you've provided
but use the |
type |
Force interpretation of the type of expression given
as |
gradient |
Control gradient derivation. If |
fixed |
An optional list of values that can be used within the DSL code. Anything you provide here is available for your calculations. In the interest of future compatibility, we check currently that all elements are scalars. In future this may become more flexible and allow passing environments, etc. Once provided, these values cannot be changed without rebuilding the model; they are fixed data. You might use these for hyperparameters that are fixed across a set of model runs, for example. |
A monty_model object derived from the expressions you provide.
# 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)")# 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)")
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.
monty_dsl_distributions()monty_dsl_distributions()
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.
monty_dsl_distributions()monty_dsl_distributions()
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.
monty_dsl_error_explain(code, how = "pretty")monty_dsl_error_explain(code, how = "pretty")
code |
The error code, as a string, in the form |
how |
How to explain the error. Options are |
Nothing, this is called for its side effect only
monty_dsl_error_explain("E201")monty_dsl_error_explain("E201")
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.
monty_dsl_parse_distribution(expr, name = NULL)monty_dsl_parse_distribution(expr, name = NULL)
expr |
An expression |
name |
Name for the expression, used in constructing messages that you can use in errors. |
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 monty_rng object (see monty_rng_create)
# A successful match monty_dsl_parse_distribution(quote(Normal(0, 1))) # An unsuccessful match monty_dsl_parse_distribution(quote(Normal()))# A successful match monty_dsl_parse_distribution(quote(Normal(0, 1))) # An unsuccessful match monty_dsl_parse_distribution(quote(Normal()))
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.
monty_example(name, ...)monty_example(name, ...)
name |
Name of the example, as a string. See Details for supported models. |
... |
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. |
A monty_model object
banannaThe 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
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.
monty_example("banana") monty_example("gaussian", diag(2))monty_example("banana") monty_example("gaussian", diag(2))
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.
monty_model(model, properties = NULL)monty_model(model, properties = NULL)
model |
A list or environment with elements as described in Details. |
properties |
Optionally, a monty_model_properties object, used to enforce or clarify properties of the model. |
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 (see
monty_rng_create). 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 a
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).
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
monty_model_function, which provides a simple interface for creating models from functions.
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.
monty_model_combine(a, b, properties = NULL, name_a = "a", name_b = "b")monty_model_combine(a, b, properties = NULL, name_a = "a", name_b = "b")
a |
The first model |
b |
The second model |
properties |
A monty_model_properties object, used to control (or enforce) properties of the combined model. |
name_a |
Name of the first model (defaulting to 'a'); you can use this to make error messages nicer to read, but it has no other practical effect. |
name_b |
Name of the first model (defaulting to 'b'); you can use this to make error messages nicer to read, but it has no other practical effect. |
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.
A monty_model object
# 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))# 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))
Compute log density for a model. This is a wrapper around the
$density property within a monty_model object.
monty_model_density(model, parameters)monty_model_density(model, parameters)
model |
A monty_model object |
parameters |
A vector or matrix of parameters |
A log-density value, or vector of log-density values
monty_model_gradient for computing gradients and monty_model_direct_sample for sampling from a model.
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))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))
Directly sample from a model. Not all models support this, and an error will be thrown if it is not possible.
monty_model_direct_sample(model, rng, named = FALSE)monty_model_direct_sample(model, rng, named = FALSE)
model |
A monty_model object |
rng |
Random number state, created by monty_rng_create.
Use of an RNG with more than one stream may or may not work as
expected; this is something we need to tidy up ( |
named |
Logical, indicating if the output should be named using the parameter names. |
A vector or matrix of sampled parameters
m <- monty_example("banana") r <- monty_rng_create() monty_model_direct_sample(m, r) monty_model_direct_sample(m, r, named = TRUE) r <- monty_rng_create(n_streams = 3) monty_model_direct_sample(m, r) monty_model_direct_sample(m, r, named = TRUE)m <- monty_example("banana") r <- monty_rng_create() monty_model_direct_sample(m, r) monty_model_direct_sample(m, r, named = TRUE) r <- monty_rng_create(n_streams = 3) monty_model_direct_sample(m, r) monty_model_direct_sample(m, r, named = TRUE)
monty_model from a function computing densityCreate 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.
monty_model_function( density, packer = NULL, fixed = NULL, allow_multiple_parameters = FALSE )monty_model_function( density, packer = NULL, fixed = NULL, allow_multiple_parameters = FALSE )
density |
A function to compute log density. It can take any number of parameters |
packer |
Optionally, a monty_packer object to control how your function parameters are packed into a numeric vector. You can typically omit this if all the arguments to your functions are present in your numeric vector and if they are all scalars. |
fixed |
Optionally, a named list of fixed values to
substitute into the call to |
allow_multiple_parameters |
Logical, indicating if passing in
vectors for all parameters will return a vector of densities.
This is |
This interface will expand in future versions of monty to support gradients, stochastic models, parameter groups and simultaneous calculation of density.
A monty_model object that computes log density with the
provided density function, given a numeric vector argument
representing all parameters.
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))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))
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.
monty_model_gradient(model, parameters, named = FALSE)monty_model_gradient(model, parameters, named = FALSE)
model |
A monty_model object |
parameters |
A vector or matrix of parameters |
named |
Logical, indicating if the output should be named using the parameter names. |
A vector or matrix of gradients
monty_model_density for log density, and monty_model_direct_sample to sample from a model
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))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))
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.
monty_model_properties( has_gradient = NULL, has_direct_sample = NULL, is_stochastic = NULL, has_parameter_groups = NULL, has_observer = NULL, allow_multiple_parameters = NULL )monty_model_properties( has_gradient = NULL, has_direct_sample = NULL, is_stochastic = NULL, has_parameter_groups = NULL, has_observer = NULL, allow_multiple_parameters = NULL )
has_gradient |
Logical, indicating if the model has a
|
has_direct_sample |
Logical, indicating if the model has a
|
is_stochastic |
Logical, indicating if the model is
stochastic. Stochastic models must supply |
has_parameter_groups |
Logical, indicating that the model can
be decomposed into parameter groups which are independent of
each other. This is indicated by using the |
has_observer |
Logical, indicating if the model has an
"observation" function, which we will describe more fully soon.
An observer is a function |
allow_multiple_parameters |
Logical, indicating if the
density calculation can support being passed a matrix of
parameters (with each column corresponding to a different
parameter set) and return a vector of densities. If |
A list of class monty_model_properties which should
not be modified.
# Default properties: monty_model_properties() # Set some properties: monty_model_properties(has_gradient = TRUE, is_stochastic = FALSE)# Default properties: monty_model_properties() # Set some properties: monty_model_properties(has_gradient = TRUE, is_stochastic = FALSE)
Split a model that has been combined by monty_model_combine() into
its constituent parts.
monty_model_split(model, prior_first = FALSE)monty_model_split(model, prior_first = FALSE)
model |
A combined model |
prior_first |
Logical, indicating that we should require that
the model component that could be the prior is listed first. If
|
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.
An unnamed list of length 2, being the component models. If one model might be the prior it will be listed first.
Create an observer to extract additional details from your model during the sampling process.
monty_observer(observe, finalise = NULL, combine = NULL, append = NULL)monty_observer(observe, finalise = NULL, combine = NULL, append = NULL)
observe |
A function that will run with arguments |
finalise |
A function that runs after a single chain has run, and you use to simplify across all samples drawn from that chain. Takes a single argument which is the list with one set of observations per sample. |
combine |
A function that runs after all chains have run, and you use to simplify across chains. Takes a single argument, which is the list with one set of observations per chain. |
append |
A function that runs after a continuation of chain has run (via monty_sample_continue. Takes two arguments representing the fully simplified observations from the first and second chains. |
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.
An object with class monty_observer which can be
passed in to monty_sample.
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".
monty_packer(scalar = NULL, array = NULL, fixed = NULL, process = NULL)monty_packer(scalar = NULL, array = NULL, fixed = NULL, process = NULL)
scalar |
Names of scalars. This is similar for listing
elements in |
array |
A list, where names correspond to the names of arrays
and values correspond to their lengths. Multiple dimensions are
allowed (so if you provide an element with two entries these
represent dimensions of a matrix). Zero-length integer vectors
or |
fixed |
A named list of fixed data to be inserted into the final unpacked list; these will be added into the final list directly. In the parameter packer context, these typically represent additional pieces of data that your model needs to run, but which you are not performing inference on. |
process |
An arbitrary R function that will be passed the
final assembled list; it may create any additional entries,
which will be concatenated onto the original list. If you use
this you should take care not to return any values with the same
names as entries listed in |
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.
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!
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.
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.
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.
# 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))# 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))
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.
monty_packer_grouped( groups, scalar = NULL, array = NULL, fixed = NULL, process = NULL, shared = NULL )monty_packer_grouped( groups, scalar = NULL, array = NULL, fixed = NULL, process = NULL, shared = NULL )
groups |
A character vector of group names. These must not
be present within any of your |
scalar |
Names of scalars. This is similar for listing
elements in |
array |
A list, where names correspond to the names of arrays
and values correspond to their lengths. Multiple dimensions are
allowed (so if you provide an element with two entries these
represent dimensions of a matrix). Zero-length integer vectors
or |
fixed |
A named list of fixed data to be inserted into the final unpacked list; these will be added into the final list directly. In the parameter packer context, these typically represent additional pieces of data that your model needs to run, but which you are not performing inference on. |
process |
An arbitrary R function that will be passed the
final assembled list for each group; it may create any
additional entries, which will be concatenated onto the
original list. If you use this you should take care not to
return any values with the same names as entries listed in
|
shared |
Names of the elements in |
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 |
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.
p <- monty_packer_grouped(c("x", "y"), c("a", "b", "c", "d", "e"), shared = c("b", "c")) p$names() p$unpack(1:8)p <- monty_packer_grouped(c("x", "y"), c("a", "b", "c", "d", "e"), shared = c("b", "c")) p$names() p$unpack(1:8)
Sample from the beta distribution
monty_random_beta(a, b, state) monty_random_n_beta(n_samples, a, b, state)monty_random_beta(a, b, state) monty_random_n_beta(n_samples, a, b, state)
a, b
|
The shape parameters |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Sample from a beta-binomial distribution. There are two parameterisations available one in terms of probability and dispersion and the other in terms of two shape parameters.
monty_random_beta_binomial_prob(size, prob, rho, state) monty_random_n_beta_binomial_prob(n_samples, size, prob, rho, state) monty_random_beta_binomial_ab(size, a, b, state) monty_random_n_beta_binomial_ab(n_samples, size, a, b, state)monty_random_beta_binomial_prob(size, prob, rho, state) monty_random_n_beta_binomial_prob(n_samples, size, prob, rho, state) monty_random_beta_binomial_ab(size, a, b, state) monty_random_n_beta_binomial_ab(n_samples, size, a, b, state)
size |
The number of trials (zero or more) |
prob |
The mean probability of success on each trial (between 0 and 1) |
rho |
The dispersion parameter (between 0 and 1) |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
a |
The first shape parameter (zero or more) |
b |
The second shape parameter (zero or more) |
A vector of random numbers, the same length as the number
of streams in state.
Sample from the binomial distribution
monty_random_binomial(size, prob, state) monty_random_n_binomial(n_samples, size, prob, state)monty_random_binomial(size, prob, state) monty_random_n_binomial(n_samples, size, prob, state)
size |
The number of trials |
prob |
The probability of success on each trial |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
state <- monty_rng_create() monty_random_binomial(10, 0.3, state) table(monty_random_n_binomial(2000, 10, 0.3, state))state <- monty_rng_create() monty_random_binomial(10, 0.3, state) table(monty_random_n_binomial(2000, 10, 0.3, state))
Sample from the Cauchy distribution
monty_random_cauchy(location, scale, state) monty_random_n_cauchy(n_samples, location, scale, state)monty_random_cauchy(location, scale, state) monty_random_n_cauchy(n_samples, location, scale, state)
location |
Location of the distribution (the same as the median and mode) |
scale |
A scale parameter which specifies the half-width at half-maximum (HWHM) |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Sample from an exponential distribution. There are two parameterisations here, one in terms of the rate of the exponential, and one in terms of the mean (or scale).
monty_random_exponential_rate(rate, state) monty_random_n_exponential_rate(n_samples, rate, state) monty_random_exponential_mean(mean, state) monty_random_n_exponential_mean(n_samples, mean, state)monty_random_exponential_rate(rate, state) monty_random_n_exponential_rate(n_samples, rate, state) monty_random_exponential_mean(mean, state) monty_random_n_exponential_mean(n_samples, mean, state)
rate |
The rate of the exponential |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
mean |
The mean of the exponential distribution (i.e., |
A vector of random numbers, the same length as the number
of streams in state.
state <- monty_rng_create() monty_random_exponential_rate(0.2, state) summary(monty_random_n_exponential_rate(2000, 0.2, state))state <- monty_rng_create() monty_random_exponential_rate(0.2, state) summary(monty_random_n_exponential_rate(2000, 0.2, state))
Sample from gamma distribution
monty_random_gamma_scale(shape, scale, state) monty_random_n_gamma_scale(n_samples, shape, scale, state) monty_random_gamma_rate(shape, rate, state) monty_random_n_gamma_rate(n_samples, shape, rate, state)monty_random_gamma_scale(shape, scale, state) monty_random_n_gamma_scale(n_samples, shape, scale, state) monty_random_gamma_rate(shape, rate, state) monty_random_n_gamma_rate(n_samples, shape, rate, state)
shape |
Shape |
scale |
Scale ' |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
rate |
Rate |
A vector of random numbers, the same length as the number
of streams in state.
Sample from a hypergeometric distribution.
monty_random_hypergeometric(n1, n2, k, state) monty_random_n_hypergeometric(n_samples, n1, n2, k, state)monty_random_hypergeometric(n1, n2, k, state) monty_random_n_hypergeometric(n_samples, n1, n2, k, state)
n1 |
The number of white balls in the urn (called n in R's rhyper) |
n2 |
The number of black balls in the urn (called m in R's rhyper) |
k |
The number of balls to draw |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Sample from a log-normal distribution
monty_random_log_normal(meanlog, sdlog, state) monty_random_n_log_normal(n_samples, meanlog, sdlog, state)monty_random_log_normal(meanlog, sdlog, state) monty_random_n_log_normal(n_samples, meanlog, sdlog, state)
meanlog |
The mean of the distribution on the log scale |
sdlog |
The standard deviation of the distribution on the log scale |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Sample from a negative binomial distribution
monty_random_negative_binomial_prob(size, prob, state) monty_random_n_negative_binomial_prob(n_samples, size, prob, state) monty_random_negative_binomial_mu(size, mu, state) monty_random_n_negative_binomial_mu(n_samples, size, mu, state)monty_random_negative_binomial_prob(size, prob, state) monty_random_n_negative_binomial_prob(n_samples, size, prob, state) monty_random_negative_binomial_mu(size, mu, state) monty_random_n_negative_binomial_mu(n_samples, size, mu, state)
size |
The target number of successful trials (zero or more) |
prob |
The probability of success on each trial (between 0 and 1) |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
mu |
The mean (zero or more) |
A vector of random numbers, the same length as the number
of streams in state.
Sample from a normal distribution
monty_random_normal(mean, sd, state, algorithm = "box_muller") monty_random_n_normal(n_samples, mean, sd, state, algorithm = "box_muller")monty_random_normal(mean, sd, state, algorithm = "box_muller") monty_random_n_normal(n_samples, mean, sd, state, algorithm = "box_muller")
mean |
The mean of the normal distribution |
sd |
The standard deviation of the normal distribution |
state |
The random number state, from monty_rng_create |
algorithm |
The algorithm to use for the normal samples;
currently |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Sample from the Poisson distribution
monty_random_poisson(lambda, state) monty_random_n_poisson(n_samples, lambda, state)monty_random_poisson(lambda, state) monty_random_n_poisson(n_samples, lambda, state)
lambda |
The mean (zero or more, length 1 or n). Only valid for lambda <= 10^7 |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Generate a random number uniformly sampled on the range 0 to 1; this is the most basic of all random number functions in monty and all other algorithms are composed from this. Quite often, you will want a number on $[0, 1]$ (e.g., for a Bernoulli trial), and this function is the most efficient way of generating one.
monty_random_real(state) monty_random_n_real(n_samples, state)monty_random_real(state) monty_random_n_real(n_samples, state)
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
state <- monty_rng_create() monty_random_real(state) monty_random_n_real(5, state)state <- monty_rng_create() monty_random_real(state) monty_random_n_real(5, state)
Sample from a truncated normal distribution
monty_random_truncated_normal(mean, sd, min, max, state) monty_random_n_truncated_normal(n_samples, mean, sd, min, max, state)monty_random_truncated_normal(mean, sd, min, max, state) monty_random_n_truncated_normal(n_samples, mean, sd, min, max, state)
mean |
The mean of the parent (untruncated) normal
distribution (this is not necessarily the mean of the truncated
distribution, unless the |
sd |
The standard deviation of the parent distribution (this is not the same as the standard deviation of the truncated distribution with finite bounds). |
min |
The lower bound (can be |
max |
The upper bound (can be |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Sample from a uniform distribution
monty_random_uniform(min, max, state) monty_random_n_uniform(n_samples, min, max, state)monty_random_uniform(min, max, state) monty_random_n_uniform(n_samples, min, max, state)
min |
The minimum value of the uniform distribution |
max |
The maximum value of the uniform distribution |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Sample from a Weibull distribution
monty_random_weibull(shape, scale, state) monty_random_n_weibull(n_samples, shape, scale, state)monty_random_weibull(shape, scale, state) monty_random_n_weibull(n_samples, shape, scale, state)
shape |
Shape |
scale |
Scale |
state |
The random number state, from monty_rng_create |
n_samples |
The number of samples to take, per stream. When using the multiple-sample interface, all other parameters are held constant (per stream). |
A vector of random numbers, the same length as the number
of streams in state.
Create a monty random number generator. This allows you to sample
random numbers from the same random number algorithms as monty
provides via C++ to dust2, and
which it uses within its samplers and filters. This function
creates the internal state, and will be passed to actual
generation functions monty_random_*, such as
monty_random_real().
monty_rng_create( n_streams = 1L, seed = NULL, n_threads = 1L, deterministic = FALSE, preserve_stream_dimension = FALSE )monty_rng_create( n_streams = 1L, seed = NULL, n_threads = 1L, deterministic = FALSE, preserve_stream_dimension = FALSE )
n_streams |
The number of streams to create (see Details) |
seed |
The initial seed of the random number generator. This
can be |
n_threads |
The number of threads to use, if OpenMP is enabled. |
deterministic |
Logical, indicating if we should use "deterministic" mode where distributions return their expectations and the state is never changed. |
preserve_stream_dimension |
Logical, indicating if the stream
dimension should be preserved in the case where |
Monty's random number generation is very different to that in base R. We have the concept of "streams" of random numbers, with a generator having 1 or many streams. Each stream is statistically independent, and can be sampled from simultaneously. When you use any of the random number functions from R, you will draw one number per stream.
Because the random number state can have multiple streams, and we
return a vector, we have a separate set of functions where
multiple numbers are requested per stream; these are all
prefixed monty_random_n_ (e.g., monty_random_n_real()). These
will return a matrix where you have used multiple streams, with
each column representing a stream. If you have a single stream
and you set preserve_stream_dimension to FALSE then we will
drop this dimension and return a matrix.
An object of class monty_rng_state, which can be
passed as the state argument to random-number producing
functions, such as monty_random_real
state <- monty_rng_create() state monty_random_real(state)state <- monty_rng_create() state monty_random_real(state)
Jump random number state. There are two "lengths" of jumps; a normal jump and a long jump. The normal jump is the distance between streams within a random number state, so if you have a multi-stream rng this shifts states left. The long jump is used to create distributed states. We will properly explain all this once the interface stabilises.
monty_rng_jump(state, n = 1) monty_rng_long_jump(state, n = 1)monty_rng_jump(state, n = 1) monty_rng_long_jump(state, n = 1)
state |
Either a |
n |
The number of jumps to take (integer, 1 or more) |
The monty_rng_state object (invisibly, modified in
place) or a raw vector, matching the input argument state
(visibly).
Get and set internal random number state
monty_rng_state(state) monty_rng_set_state(value, state)monty_rng_state(state) monty_rng_set_state(value, state)
state |
The random number state, from monty_rng_create |
value |
A vector of raw values, typically the result of
exporting a random state with |
A vector of raws
s1 <- monty_rng_create() r1 <- monty_rng_state(s1) s2 <- monty_rng_create(seed = r1) identical(r1, monty_rng_state(s2)) monty_random_real(s1) monty_random_real(s2) monty_rng_set_state(r1, s1) monty_random_real(s1) monty_random_real(s1) monty_random_real(s2)s1 <- monty_rng_create() r1 <- monty_rng_state(s1) s2 <- monty_rng_create(seed = r1) identical(r1, monty_rng_state(s2)) monty_random_real(s1) monty_random_real(s2) monty_rng_set_state(r1, s1) monty_random_real(s1) monty_random_real(s1) monty_random_real(s2)
callr
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.
monty_runner_callr(n_workers, progress = NULL)monty_runner_callr(n_workers, progress = NULL)
n_workers |
The number of workers to use. This should be no
larger than the number of chains (though this is harmless) and
no larger than the total number of cores available on your
computer. Ideally the number of chains you want to run is a
multiple of this number (for example, if you had 8 chains, then
1, 2, 4, and 8 are good choices of |
progress |
Optional logical, indicating if we should print a
progress bar while running. If |
A runner of class monty_runner that can be passed to
monty_sample()
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.
monty_runner_parallel(n_workers)monty_runner_parallel(n_workers)
n_workers |
Number of workers to create a cluster from. In a
multi-user setting be careful not to set this to more cores than
you are allowed to use. You can use |
A runner of class monty_runner that can be passed to
monty_sample()
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).
monty_runner_serial(progress = NULL)monty_runner_serial(progress = NULL)
progress |
Optional logical, indicating if we should print a
progress bar while running. If |
A runner of class monty_runner that can be passed to
monty_sample()
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.
monty_runner_simultaneous(progress = NULL)monty_runner_simultaneous(progress = NULL)
progress |
Optional logical, indicating if we should print a
progress bar while running. If |
A runner of class monty_runner that can be passed to
monty_sample()
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)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)
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.
monty_sample( model, sampler, n_steps, initial = NULL, n_chains = 1L, runner = NULL, restartable = FALSE, burnin = NULL, thinning_factor = NULL )monty_sample( model, sampler, n_steps, initial = NULL, n_chains = 1L, runner = NULL, restartable = FALSE, burnin = NULL, thinning_factor = NULL )
model |
The model to sample from; this should be a
|
sampler |
A sampler to use. These will be described later,
but we hope to make these reasonably easy to implement so that
we can try out different sampling ideas. For now, the only
sampler implemented is |
n_steps |
The number of steps to run the sampler for. |
initial |
Optionally, initial parameter values for the
sampling. If not given, we sample from the model (or its
prior). Alternatively, you can provide a |
n_chains |
Number of chains to run. The default is to run a single chain, but you will likely want to run more. |
runner |
A runner for your chains. The default option is to
run chains in series (via monty_runner_serial). The only
other current option is monty_runner_parallel which uses the
|
restartable |
Logical, indicating if the chains should be restartable. This will add additional data to the chains object. |
burnin |
Number of steps to discard as burnin. This affects
only the recording of steps as your chains run; we don't record
the first |
thinning_factor |
A thinning factor to apply while the chain
is running. If given, then we save every |
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.
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.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.
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.
monty_sample_continue( samples, n_steps, restartable = FALSE, runner = NULL, append = TRUE )monty_sample_continue( samples, n_steps, restartable = FALSE, runner = NULL, append = TRUE )
samples |
A |
n_steps |
The number of new steps to run |
restartable |
Logical, indicating if the chains should be restartable. This will add additional data to the chains object. |
runner |
Optionally, a runner for your chains. The default
is to continue with the backend that you used to start the
chains via monty_sample (or on the previous restart with
this function). You can use this argument to change the runner,
which might be useful if transferring a pilot run from a
high-resource environment to a lower-resource environment. If
given, must be a |
append |
Logical, indicating if we should append the results of the resumed chain together with the original chain. |
A list of parameters and densities
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.
monty_sample_manual_cleanup(path)monty_sample_manual_cleanup(path)
path |
The path used in the call to monty_sample_manual_prepare |
Nothing, called for side effects only.
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)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)
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.
monty_sample_manual_collect( path, samples = NULL, restartable = FALSE, append = TRUE )monty_sample_manual_collect( path, samples = NULL, restartable = FALSE, append = TRUE )
path |
The path used in the call to monty_sample_manual_prepare |
samples |
Samples from the parent run. You need to provide
these where |
restartable |
Logical, indicating if the chains should be restartable. This will add additional data to the chains object. Note that this is controlled at chain collection and not creation. |
append |
Logical, indicating if we should append the results of the resumed chain together with the original chain. |
A monty_samples object.
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)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)
Get information about the status of manually scheduled samples.
monty_sample_manual_info(path)monty_sample_manual_info(path)
path |
The path used in the call to monty_sample_manual_prepare |
Invisibly, a logical vector, TRUE for completed chains
and FALSE for incomplete chains.
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)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)
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.
monty_sample_manual_prepare( model, sampler, n_steps, path, initial = NULL, n_chains = 1L, burnin = NULL, thinning_factor = NULL )monty_sample_manual_prepare( model, sampler, n_steps, path, initial = NULL, n_chains = 1L, burnin = NULL, thinning_factor = NULL )
model |
The model to sample from; this should be a
|
sampler |
A sampler to use. These will be described later,
but we hope to make these reasonably easy to implement so that
we can try out different sampling ideas. For now, the only
sampler implemented is |
n_steps |
The number of steps to run the sampler for. |
path |
The path to write inputs and outputs to. This should
be a path to a directory which does not yet exist, or which is
empty; we will create one here. The contents of this directory
are managed by |
initial |
Optionally, initial parameter values for the
sampling. If not given, we sample from the model (or its
prior). Alternatively, you can provide a |
n_chains |
Number of chains to run. The default is to run a single chain, but you will likely want to run more. |
burnin |
Number of steps to discard as burnin. This affects
only the recording of steps as your chains run; we don't record
the first |
thinning_factor |
A thinning factor to apply while the chain
is running. If given, then we save every |
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.
Invisibly, the path used to store files (the same as the
value of the path argument)
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.
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)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)
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.
monty_sample_manual_prepare_continue( samples, n_steps, path, save_samples = "hash" )monty_sample_manual_prepare_continue( samples, n_steps, path, save_samples = "hash" )
samples |
A |
n_steps |
The number of steps to run the sampler for. |
path |
The path to write inputs and outputs to. This should
be a path to a directory which does not yet exist, or which is
empty; we will create one here. The contents of this directory
are managed by |
save_samples |
Control over saving samples into the inputs.
The choices here are |
Invisibly, the path used to store files (the same as the
value of the path argument)
Run a chain that was prepared using monty_sample_manual_prepare.
monty_sample_manual_run(chain_id, path, progress = NULL)monty_sample_manual_run(chain_id, path, progress = NULL)
chain_id |
The id for the chain to run, an integer. If you
provide an integer that does not correspond to a chain in 1 to
|
path |
The path used in the call to monty_sample_manual_prepare |
progress |
Optional logical, indicating if we should print a
progress bar while running. If |
There is no lock mechanism; you can start a single chain many times. Don't do this.
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)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)
Create an adaptive Metropolis-Hastings sampler, which will tune its variance covariance matrix (vs the simple random walk sampler monty_sampler_random_walk).
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" )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" )
initial_vcv |
An initial variance covariance matrix; we'll start using this in the proposal, which will gradually become more weighted towards the empirical covariance matrix calculated from the chain. |
initial_vcv_weight |
Weight of the initial variance-covariance matrix used to build the proposal of the random-walk. Higher values translate into higher confidence of the initial variance-covariance matrix and means that update from additional samples will be slower. |
initial_scaling |
The initial scaling of the variance covariance matrix to be used to generate the multivariate normal proposal for the random-walk Metropolis-Hastings algorithm. To generate the proposal matrix, the weighted variance covariance matrix is multiplied by the scaling parameter squared times 2.38^2 / n_pars (where n_pars is the number of fitted parameters). Thus, in a Gaussian target parameter space, the optimal scaling will be around 1. |
initial_scaling_weight |
The initial weight used in the scaling update.
The scaling weight will increase after the first |
min_scaling |
The minimum scaling of the variance covariance matrix to be used to generate the multivariate normal proposal for the random-walk Metropolis-Hastings algorithm. |
scaling_increment |
The scaling increment which is added or
subtracted to the scaling factor of the variance-covariance
after each adaptive step. If |
log_scaling_update |
Logical, whether or not changes to the scaling parameter are made on the log-scale. |
acceptance_target |
The target for the fraction of proposals that should be accepted (optimally) for the adaptive part of the chain. |
forget_rate |
The rate of forgetting early parameter sets from the
empirical variance-covariance matrix in the MCMC chains. For example,
|
forget_end |
The final iteration at which early parameter sets can
be forgotten. Setting |
adapt_end |
The final iteration at which we can adapt the multivariate normal proposal. Thereafter the empirical variance-covariance matrix, its scaling and its weight remain fixed. This allows the adaptation to be switched off at a certain point to help ensure convergence of the chain. |
pre_diminish |
The number of updates before adaptation of the scaling
parameter starts to diminish. Setting |
boundaries |
Control the behaviour of proposals that are outside the model domain. The supported options are:
The initial point selected will lie within the domain, as this is enforced by monty_sample. |
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.
A monty_sampler object, which can be used with
monty_sample
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))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))
Create a Hamiltonian Monte Carlo sampler, implemented using the leapfrog algorithm.
monty_sampler_hmc( epsilon = 0.015, n_integration_steps = 10, vcv = NULL, debug = FALSE )monty_sampler_hmc( epsilon = 0.015, n_integration_steps = 10, vcv = NULL, debug = FALSE )
epsilon |
The step size of the HMC steps |
n_integration_steps |
The number of HMC steps per step |
vcv |
A variance-covariance matrix for the momentum vector. The default uses an identity matrix. |
debug |
Logical, indicating if we should save all intermediate points and their gradients. This will add a vector "history" to the details after the integration. This will slow things down though as we accumulate the history inefficiently. |
A monty_sampler object, which can be used with
monty_sample
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.
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" )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" )
initial_vcv |
An initial variance covariance matrix; we'll start using this in the proposal, which will gradually become more weighted towards the empirical covariance matrix calculated from the chain. |
initial_vcv_weight |
Weight of the initial variance-covariance matrix used to build the proposal of the random-walk. Higher values translate into higher confidence of the initial variance-covariance matrix and means that update from additional samples will be slower. |
initial_scaling |
The initial scaling of the variance covariance matrix to be used to generate the multivariate normal proposal for the random-walk Metropolis-Hastings algorithm. To generate the proposal matrix, the weighted variance covariance matrix is multiplied by the scaling parameter squared times 2.38^2 / n_pars (where n_pars is the number of fitted parameters). Thus, in a Gaussian target parameter space, the optimal scaling will be around 1. |
initial_scaling_weight |
The initial weight used in the scaling update.
The scaling weight will increase after the first |
min_scaling |
The minimum scaling of the variance covariance matrix to be used to generate the multivariate normal proposal for the random-walk Metropolis-Hastings algorithm. |
scaling_increment |
The scaling increment which is added or
subtracted to the scaling factor of the variance-covariance
after each adaptive step. If |
log_scaling_update |
Logical, whether or not changes to the scaling parameter are made on the log-scale. |
acceptance_target |
The target for the fraction of proposals that should be accepted (optimally) for the adaptive part of the chain. |
forget_rate |
The rate of forgetting early parameter sets from the
empirical variance-covariance matrix in the MCMC chains. For example,
|
forget_end |
The final iteration at which early parameter sets can
be forgotten. Setting |
adapt_end |
The final iteration at which we can adapt the multivariate normal proposal. Thereafter the empirical variance-covariance matrix, its scaling and its weight remain fixed. This allows the adaptation to be switched off at a certain point to help ensure convergence of the chain. |
pre_diminish |
The number of updates before adaptation of the scaling
parameter starts to diminish. Setting |
boundaries |
Control the behaviour of proposals that are outside the model domain. The supported options are:
The initial point selected will lie within the domain, as this is enforced by monty_sample. |
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.
A monty_sampler object, which can be used with
monty_sample
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.
monty_sampler_nested_random_walk(vcv, boundaries = "reflect")monty_sampler_nested_random_walk(vcv, boundaries = "reflect")
vcv |
A list of variance covariance matrices. We expect this
to be a list with elements |
boundaries |
Control the behaviour of proposals that are outside the model domain. The supported options are:
The initial point selected will lie within the domain, as this is enforced by monty_sample. |
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.
A monty_sampler object, which can be used with
monty_sample
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).
monty_sampler_random_walk( vcv = NULL, boundaries = "reflect", rerun_every = Inf, rerun_random = TRUE )monty_sampler_random_walk( vcv = NULL, boundaries = "reflect", rerun_every = Inf, rerun_random = TRUE )
vcv |
A variance covariance matrix for the proposal. |
boundaries |
Control the behaviour of proposals that are outside the model domain. The supported options are:
The initial point selected will lie within the domain, as this is enforced by monty_sample. |
rerun_every |
Optional integer giving the frequency at which
we should rerun the model on the current "accepted" parameters to
obtain a new density for stochastic models. The default for this
( |
rerun_random |
Logical, controlling the behaviour of
rerunning (when |
A monty_sampler object, which can be used with
monty_sample
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.
monty_samples_thin(samples, thinning_factor = NULL, burnin = NULL)monty_samples_thin(samples, thinning_factor = NULL, burnin = NULL)
samples |
A |
thinning_factor |
Optional integer thinning factor. If given,
then we save every |
burnin |
Number of steps to discard as burnin from the start of the chain. |
A monty_samples object (as for monty_sample()),
typically with fewer samples.
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.
Trace calls to R's random-number-generating functions, to detect unexpected use of random number generation outside of monty's control.
with_trace_random(code, max_calls = 5, show_stack = FALSE)with_trace_random(code, max_calls = 5, show_stack = FALSE)
code |
Code to run with tracing on |
max_calls |
Maximum number of calls to report. The default is 5 |
show_stack |
Logical, indicating if we should show the stack at the point of the call |
The result of evaluating code