diff --git a/DESCRIPTION b/DESCRIPTION index 63518ca8..528079fe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.2.19 +Version: 0.2.22 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), @@ -31,7 +31,6 @@ Suggests: coda, decor, knitr, - dust, mockery, mvtnorm, numDeriv, @@ -43,5 +42,3 @@ Suggests: Config/testthat/edition: 3 Language: en-GB VignetteBuilder: knitr -Remotes: - mrc-ide/dust diff --git a/R/cpp11.R b/R/cpp11.R index 9c38a8cf..15a17c9c 100644 --- a/R/cpp11.R +++ b/R/cpp11.R @@ -40,8 +40,12 @@ monty_rng_binomial <- function(ptr, n, r_size, r_prob, n_threads, is_float) { .Call(`_monty_monty_rng_binomial`, ptr, n, r_size, r_prob, n_threads, is_float) } -monty_rng_nbinomial <- function(ptr, n, r_size, r_prob, n_threads, is_float) { - .Call(`_monty_monty_rng_nbinomial`, ptr, n, r_size, r_prob, n_threads, is_float) +monty_rng_negative_binomial_prob <- function(ptr, n, r_size, r_prob, n_threads, is_float) { + .Call(`_monty_monty_rng_negative_binomial_prob`, ptr, n, r_size, r_prob, n_threads, is_float) +} + +monty_rng_negative_binomial_mu <- function(ptr, n, r_size, r_mu, n_threads, is_float) { + .Call(`_monty_monty_rng_negative_binomial_mu`, ptr, n, r_size, r_mu, n_threads, is_float) } monty_rng_hypergeometric <- function(ptr, n, r_n1, r_n2, r_k, n_threads, is_float) { @@ -76,6 +80,18 @@ monty_rng_state <- function(ptr, is_float) { .Call(`_monty_monty_rng_state`, ptr, is_float) } +cpp_monty_random_real <- function(ptr) { + .Call(`_monty_cpp_monty_random_real`, ptr) +} + +cpp_monty_random_binomial <- function(r_size, r_prob, ptr) { + .Call(`_monty_cpp_monty_random_binomial`, r_size, r_prob, ptr) +} + +cpp_monty_random_exponential_rate <- function(r_rate, ptr) { + .Call(`_monty_cpp_monty_random_exponential_rate`, r_rate, ptr) +} + monty_rng_pointer_init <- function(n_streams, seed, long_jump, algorithm) { .Call(`_monty_monty_rng_pointer_init`, n_streams, seed, long_jump, algorithm) } diff --git a/R/domain.R b/R/domain.R index 6c963886..a33a7477 100644 --- a/R/domain.R +++ b/R/domain.R @@ -54,7 +54,7 @@ monty_domain_expand <- function(domain, packer) { arg = "domain") } - nms_full <- packer$parameters + nms_full <- packer$names() nms_map <- packer$unpack(nms_full) nms_logical <- names(nms_map) @@ -76,5 +76,5 @@ monty_domain_expand <- function(domain, packer) { domain[!i, , drop = FALSE]) } - domain[order(match(rownames(domain), packer$parameters)), , drop = FALSE] + domain[order(match(rownames(domain), packer$names())), , drop = FALSE] } diff --git a/R/dsl-distributions.R b/R/dsl-distributions.R index d8af6aae..6b9ee832 100644 --- a/R/dsl-distributions.R +++ b/R/dsl-distributions.R @@ -103,7 +103,7 @@ distr_gamma_rate <- distribution( distr_gamma_scale <- distribution( name = "Gamma", - variant = "rate", + variant = "scale", density = function(x, shape, scale) { dgamma(x, shape, scale = scale, log = TRUE) }, @@ -125,6 +125,36 @@ distr_hypergeometric <- distribution( mean = quote(k * n1 / (n1 + n2))), cpp = list(density = "hypergeometric", sample = "hypergeometric")) +distr_negative_binomial_prob <- distribution( + name = "NegativeBinomial", + variant = "prob", + density = function(x, size, prob) { + dnbinom(x, size, prob = prob, log = TRUE) + }, + domain = c(0, Inf), + sample = function(rng, size, prob) rng$negative_binomial_prob(1, size, prob), + expr = list( + density = quote(lgamma(x + size) - lgamma(size) - lgamma(x + 1) + + x * log(1 - prob) + size * log(prob)), + mean = quote(size * (1 - prob) / prob)), + cpp = list(density = "negative_binomial_prob", + sample = "negative_binomial_prob")) + +distr_negative_binomial_mu <- distribution( + name = "NegativeBinomial", + variant = "mu", + density = function(x, size, mu) { + dnbinom(x, size, mu = mu, log = TRUE) + }, + domain = c(0, Inf), + sample = function(rng, size, mu) rng$negative_binomial_mu(1, size, mu), + expr = list( + density = quote(lgamma(x + size) - lgamma(size) - lgamma(x + 1) + + size * log(size) + x * log(mu) - + (size + x) * log(size + mu)), + mean = quote(mu)), + cpp = list(density = "negative_binomial_mu", sample = "negative_binomial_mu")) + distr_normal <- distribution( name = "Normal", density = function(x, mean, sd) dnorm(x, mean, sd, log = TRUE), @@ -166,6 +196,8 @@ dsl_distributions <- local({ distr_gamma_rate, distr_gamma_scale, distr_hypergeometric, + distr_negative_binomial_prob, + distr_negative_binomial_mu, distr_normal, distr_poisson, distr_uniform) diff --git a/R/dsl-generate.R b/R/dsl-generate.R index 01ac83d7..4b8942b6 100644 --- a/R/dsl-generate.R +++ b/R/dsl-generate.R @@ -56,7 +56,7 @@ dsl_generate_direct_sample <- function(dat, env, meta) { exprs <- lapply(dat$exprs, dsl_generate_sample_expr, meta) body <- c(call("<-", meta[["pars"]], quote(list())), exprs, - bquote(unlist(.(meta[["pars"]])[packer$parameters], FALSE, FALSE))) + bquote(unlist(.(meta[["pars"]])[packer$names()], FALSE, FALSE))) as_function(alist(rng = ), body, env) } diff --git a/R/model-function.R b/R/model-function.R index af52af7e..1e2cc5c1 100644 --- a/R/model-function.R +++ b/R/model-function.R @@ -63,7 +63,7 @@ monty_model_function <- function(density, packer = NULL, fixed = NULL) { } monty_model( - list(parameters = packer$parameters, + list(parameters = packer$names(), density = function(x) { rlang::inject(density(!!!packer$unpack(x))) })) diff --git a/R/packer.R b/R/packer.R index daa7ebe8..bc84754b 100644 --- a/R/packer.R +++ b/R/packer.R @@ -1,9 +1,11 @@ -##' Build a parameter packer, which can be used in models to translate -##' between an unstructured vector of numbers (the vector being -##' updated by an MCMC for example) to a structured list of named -##' values, which is easier to program against. 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 +##' 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 models'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". ##' ##' There are several places where it is most convenient to work in an @@ -145,57 +147,64 @@ ##' vector vs packing an array where all dimensions are 1. See the ##' examples, and please let us know if the behaviour needs changing. ##' -##' @title Build a parameter packer +##' @title Build a packer ##' -##' @param scalar Names of scalar parameters. This is similar for -##' listing elements in `array` with values of 1, though elements in +##' @param scalar Names of scalars. This is similar for listing +##' elements in `array` with values of 1, though elements in ##' `scalar` will be placed ahead of those listed in `array` within ##' the final parameter vector, and elements in `array` will have ##' generated names that include square brackets. ##' -##' @param array A list, where names correspond to the names of array -##' parameters and values correspond to the lengths of parameters. -##' Multiple dimensions are allowed (so if you provide an element -##' with two entries these represent dimensions of a matrix). -##' Zero-length integer vectors or `NULL` values are counted as -##' scalars, which allows you to put scalars at positions other than -##' the front of the packing vector. In future, you may be able to -##' use *strings* as values for the lengths, in which case these -##' will be looked for within `fixed`. -##' -##' @param fixed A named list of fixed parameters; these will be added -##' into the final list directly. These typically represent -##' additional pieces of data that your model needs to run, but -##' which you are not performing inference on. +##' @param 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 `NULL` values are counted as scalars, which allows you to put +##' scalars at positions other than the front of the packing +##' vector. In future, you may be able to use *strings* as values +##' for the lengths, in which case these will be looked for within +##' `fixed`. +##' +##' @param 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. ##' ##' @param process An arbitrary R function that will be passed the -##' final assembled parameter 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 `scalar`, `array` or -##' `fixed`, as this is an error (this is so that `pack()` is -##' not broken). We will likely play around with this process in -##' future in order to get automatic differentiation to work. -##' -##' @return An object of class `monty_packer`, which has three -##' elements: -##' -##' * `parameters`: a character vector of computed parameter names; -##' these are the names that your statistical model will use. +##' 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 `scalar`, `array` or `fixed`, as this +##' is an error (this is so that `pack()` is not broken). We will +##' likely play around with this process in future in order to get +##' automatic differentiation to work. +##' +##' @return 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 -##' parameters back into a numeric vector suitable for the +##' +##' * `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. +##' `preprocess` function and present in `fixed`. +##' ##' * `index`: a function which produces a named list where each -##' element has the name of a value in `parameters` and each value -##' has the indices within an unstructured vector where these values -##' can be found. -##' * `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. +##' 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! ##' ##' @export ##' @@ -260,7 +269,7 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL, process = NULL) { call <- environment() - parameters <- character(0) + nms <- character(0) idx <- list() shape <- list() @@ -279,7 +288,7 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL, } shape[scalar] <- rep(list(integer()), length(scalar)) idx[scalar] <- as.list(seq_along(scalar)) - parameters <- c(parameters, scalar) + nms <- c(nms, scalar) } len <- length(scalar) # start arrays after scalars @@ -287,7 +296,7 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL, assert_named(array, unique = TRUE, call = call) for (nm in names(array)) { tmp <- prepare_pack_array(nm, array[[nm]], call) - parameters <- c(parameters, tmp$names) + nms <- c(nms, tmp$names) shape[[nm]] <- tmp$shape idx[[nm]] <- seq_len(tmp$n) + len len <- len + tmp$n @@ -325,9 +334,9 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL, unpack <- function(x) { if (is.null(dim(x))) { - unpack_vector(x, parameters, len, idx, shape, fixed, process) + unpack_vector(x, nms, len, idx, shape, fixed, process) } else { - unpack_array(x, parameters, len, idx, shape, fixed, process) + unpack_array(x, nms, len, idx, shape, fixed, process) } } @@ -379,7 +388,7 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL, list(index = index, packer = packer) } - ret <- list(parameters = parameters, + ret <- list(names = function() nms, unpack = unpack, pack = pack, index = function() idx, @@ -438,7 +447,7 @@ array_indices <- function(shape) { print.monty_packer <- function(x, ...) { cli::cli_h1("") cli::cli_alert_info( - "Packing {length(x$parameters)} parameter{?s}: {squote(x$parameters)}") + "Packing {length(x$names())} parameter{?s}: {squote(x$names())}") cli::cli_alert_info( "Use '$pack()' to convert from a list to a vector") cli::cli_alert_info( @@ -448,10 +457,10 @@ print.monty_packer <- function(x, ...) { } -unpack_vector <- function(x, parameters, len, idx, shape, fixed, process) { +unpack_vector <- function(x, nms, len, idx, shape, fixed, process) { call <- parent.frame() if (!is.null(names(x))) { - if (!identical(names(x), parameters)) { + if (!identical(names(x), nms)) { ## Here, we could do better I think with this message; we ## might pass thropuigh empty names, and produce some summary ## of different names. Something for later though. @@ -476,7 +485,7 @@ unpack_vector <- function(x, parameters, len, idx, shape, fixed, process) { err <- intersect(names(extra), names(res)) if (length(err) > 0) { cli::cli_abort( - c("'process()' is trying to overwrite entries in parameters", + c("'process()' is trying to overwrite entries in your list", i = paste("The 'process()' function should only create elements", "that are not already present in 'scalar', 'array'", "or 'fixed', as this lets us reverse the transformation", @@ -494,10 +503,10 @@ unpack_vector <- function(x, parameters, len, idx, shape, fixed, process) { } -unpack_array <- function(x, parameters, len, idx, shape, fixed, process) { +unpack_array <- function(x, nms, len, idx, shape, fixed, process) { call <- parent.frame() dn <- dimnames(x) - if (!is.null(dn) && !is.null(dn[[1]]) && !identical(dn[[1]], parameters)) { + if (!is.null(dn) && !is.null(dn[[1]]) && !identical(dn[[1]], nms)) { ## See comment above about reporting on this better cli::cli_abort("Incorrect rownames in input") } diff --git a/R/rng.R b/R/rng.R index 35c98919..27cf60d7 100644 --- a/R/rng.R +++ b/R/rng.R @@ -96,7 +96,10 @@ ##' rng$binomial(5, 10, 0.3) ##' ##' # Negative binomially distributed random numbers with size and prob -##' rng$nbinomial(5, 10, 0.3) +##' 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) @@ -283,9 +286,25 @@ monty_rng <- R6::R6Class( ##' (between 0 and 1, length 1 or n) ##' ##' @param n_threads Number of threads to use; see Details - nbinomial = function(n, size, prob, n_threads = 1L) { - monty_rng_nbinomial(private$ptr, n, size, prob, n_threads, - private$float) + negative_binomial_prob = function(n, size, prob, n_threads = 1L) { + monty_rng_negative_binomial_prob(private$ptr, n, size, prob, n_threads, + private$float) + }, + + ##' @description Generate `n` numbers from a negative binomial distribution + ##' + ##' @param n Number of samples to draw (per stream) + ##' + ##' @param size The target number of successful trials + ##' (zero or more, length 1 or n) + ##' + ##' @param mu The mean + ##' (zero or more, length 1 or n) + ##' + ##' @param n_threads Number of threads to use; see Details + negative_binomial_mu = function(n, size, mu, n_threads = 1L) { + monty_rng_negative_binomial_mu(private$ptr, n, size, mu, n_threads, + private$float) }, ##' @description Generate `n` numbers from a hypergeometric distribution diff --git a/R/rng2.R b/R/rng2.R new file mode 100644 index 00000000..8cded76d --- /dev/null +++ b/R/rng2.R @@ -0,0 +1,23 @@ +monty_random_alloc <- function(n_streams = 1L, seed = NULL, + deterministic = FALSE) { + float <- FALSE + ptr <- monty_rng_alloc(seed, n_streams, deterministic, float) + ret <- list(ptr = ptr, n_streams = n_streams, deterministic = deterministic) + class(ret) <- "monty_random_state" + ret +} + + +monty_random_real <- function(state) { + cpp_monty_random_real(state$ptr) +} + + +monty_random_binomial <- function(size, prob, state) { + cpp_monty_random_binomial(size, prob, state$ptr) +} + + +monty_random_exponential_rate <- function(rate, state) { + cpp_monty_random_exponential_rate(rate, state$ptr) +} diff --git a/inst/include/monty/random/nbinomial.hpp b/inst/include/monty/random/nbinomial.hpp deleted file mode 100644 index 004623fb..00000000 --- a/inst/include/monty/random/nbinomial.hpp +++ /dev/null @@ -1,41 +0,0 @@ -#pragma once - -#include - -#include "monty/random/gamma.hpp" -#include "monty/random/poisson.hpp" -#include "monty/random/generator.hpp" - -namespace monty { -namespace random { - -namespace { - -template -void nbinomial_validate(real_type size, real_type prob) { - if(!R_FINITE(size) || !R_FINITE(prob) || size <= 0 || prob <= 0 || prob > 1) { - char buffer[256]; - snprintf(buffer, 256, - "Invalid call to nbinomial with size = %g, prob = %g", - size, prob); - monty::utils::fatal_error(buffer); - } -} - -} - -template -real_type nbinomial(rng_state_type& rng_state, real_type size, real_type prob) { -#ifdef __CUDA_ARCH__ - static_assert("nbinomial() not implemented for GPU targets"); -#endif - nbinomial_validate(size, prob); - - if (rng_state.deterministic) { - return (1 - prob) * size / prob; - } - return (prob == 1) ? 0 : poisson(rng_state, gamma_scale(rng_state, size, (1 - prob) / prob)); -} - -} -} diff --git a/inst/include/monty/random/negative_binomial.hpp b/inst/include/monty/random/negative_binomial.hpp new file mode 100644 index 00000000..2c3cb61d --- /dev/null +++ b/inst/include/monty/random/negative_binomial.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include + +#include "monty/random/gamma.hpp" +#include "monty/random/poisson.hpp" +#include "monty/random/generator.hpp" + +namespace monty { +namespace random { + +namespace { + +template +void negative_binomial_validate(real_type size, real_type prob) { + if(!R_FINITE(size) || !R_FINITE(prob) || size <= 0 || prob <= 0 || prob > 1) { + char buffer[256]; + snprintf(buffer, 256, + "Invalid call to negative_binomial with size = %g, prob = %g", + size, prob); + monty::utils::fatal_error(buffer); + } +} + +} + +template +real_type negative_binomial_prob(rng_state_type& rng_state, real_type size, real_type prob) { +#ifdef __CUDA_ARCH__ + static_assert("negative_binomial_prob() not implemented for GPU targets"); +#endif + negative_binomial_validate(size, prob); + + if (rng_state.deterministic) { + return (1 - prob) * size / prob; + } + return (prob == 1) ? 0 : poisson(rng_state, gamma_scale(rng_state, size, (1 - prob) / prob)); +} + +template +real_type negative_binomial_mu(rng_state_type& rng_state, real_type size, real_type mu) { +#ifdef __CUDA_ARCH__ + static_assert("negative_binomial_mu() not implemented for GPU targets"); +#endif + const auto prob = size / (size + mu); + negative_binomial_validate(size, prob); + + return negative_binomial_prob(rng_state, size, prob); +} + +} +} diff --git a/inst/include/monty/random/random.hpp b/inst/include/monty/random/random.hpp index 508a52d8..93df6424 100644 --- a/inst/include/monty/random/random.hpp +++ b/inst/include/monty/random/random.hpp @@ -10,7 +10,7 @@ #include "monty/random/gamma.hpp" #include "monty/random/hypergeometric.hpp" #include "monty/random/multinomial.hpp" -#include "monty/random/nbinomial.hpp" +#include "monty/random/negative_binomial.hpp" #include "monty/random/normal.hpp" #include "monty/random/poisson.hpp" #include "monty/random/uniform.hpp" diff --git a/man/monty_packer.Rd b/man/monty_packer.Rd index c4ac8cc7..cd192499 100644 --- a/man/monty_packer.Rd +++ b/man/monty_packer.Rd @@ -2,70 +2,75 @@ % Please edit documentation in R/packer.R \name{monty_packer} \alias{monty_packer} -\title{Build a parameter packer} +\title{Build a packer} \usage{ monty_packer(scalar = NULL, array = NULL, fixed = NULL, process = NULL) } \arguments{ -\item{scalar}{Names of scalar parameters. This is similar for -listing elements in \code{array} with values of 1, though elements in +\item{scalar}{Names of scalars. This is similar for listing +elements in \code{array} with values of 1, though elements in \code{scalar} will be placed ahead of those listed in \code{array} within the final parameter vector, and elements in \code{array} will have generated names that include square brackets.} -\item{array}{A list, where names correspond to the names of array -parameters and values correspond to the lengths of parameters. -Multiple dimensions are allowed (so if you provide an element -with two entries these represent dimensions of a matrix). -Zero-length integer vectors or \code{NULL} values are counted as -scalars, which allows you to put scalars at positions other than -the front of the packing vector. In future, you may be able to -use \emph{strings} as values for the lengths, in which case these -will be looked for within \code{fixed}.} - -\item{fixed}{A named list of fixed parameters; these will be added -into the final list directly. These typically represent -additional pieces of data that your model needs to run, but -which you are not performing inference on.} +\item{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 \code{NULL} values are counted as scalars, which allows you to put +scalars at positions other than the front of the packing +vector. In future, you may be able to use \emph{strings} as values +for the lengths, in which case these will be looked for within +\code{fixed}.} + +\item{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.} \item{process}{An arbitrary R function that will be passed the -final assembled parameter list; it may create any \emph{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 \code{scalar}, \code{array} or -\code{fixed}, as this is an error (this is so that \code{pack()} is -not broken). We will likely play around with this process in -future in order to get automatic differentiation to work.} +final assembled list; it may create any \emph{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 \code{scalar}, \code{array} or \code{fixed}, as this +is an error (this is so that \code{pack()} is not broken). We will +likely play around with this process in future in order to get +automatic differentiation to work.} } \value{ -An object of class \code{monty_packer}, which has three -elements: +An object of class \code{monty_packer}, which has elements: \itemize{ -\item \code{parameters}: a character vector of computed parameter names; -these are the names that your statistical model will use. +\item \code{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. \item \code{unpack}: a function that can unpack an unstructured vector (say, from your statistical model parameters) into a structured list (say, for your generative model) -\item \code{pack}: a function that can pack your structured list of -parameters back into a numeric vector suitable for the +\item \code{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 -\code{preprocess} function. +\code{preprocess} function and present in \code{fixed}. \item \code{index}: a function which produces a named list where each -element has the name of a value in \code{parameters} and each value -has the indices within an unstructured vector where these values -can be found. -\item \code{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. +element has the name of a value in \code{scalar} or \code{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. +\item \code{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! } } \description{ -Build a parameter packer, which can be used in models to translate -between an unstructured vector of numbers (the vector being -updated by an MCMC for example) to a structured list of named -values, which is easier to program against. 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 +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 models'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". } \details{ diff --git a/man/monty_rng.Rd b/man/monty_rng.Rd index 6962e231..4642f365 100644 --- a/man/monty_rng.Rd +++ b/man/monty_rng.Rd @@ -97,7 +97,10 @@ rng$normal(5, 4, 2) rng$binomial(5, 10, 0.3) # Negative binomially distributed random numbers with size and prob -rng$nbinomial(5, 10, 0.3) +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) @@ -140,7 +143,8 @@ rng$multinomial(5, 10, c(0.1, 0.3, 0.5, 0.1)) \item \href{#method-monty_rng-uniform}{\code{monty_rng$uniform()}} \item \href{#method-monty_rng-normal}{\code{monty_rng$normal()}} \item \href{#method-monty_rng-binomial}{\code{monty_rng$binomial()}} -\item \href{#method-monty_rng-nbinomial}{\code{monty_rng$nbinomial()}} +\item \href{#method-monty_rng-negative_binomial_prob}{\code{monty_rng$negative_binomial_prob()}} +\item \href{#method-monty_rng-negative_binomial_mu}{\code{monty_rng$negative_binomial_mu()}} \item \href{#method-monty_rng-hypergeometric}{\code{monty_rng$hypergeometric()}} \item \href{#method-monty_rng-gamma_scale}{\code{monty_rng$gamma_scale()}} \item \href{#method-monty_rng-gamma_rate}{\code{monty_rng$gamma_rate()}} @@ -341,12 +345,12 @@ Generate \code{n} numbers from a binomial distribution } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-nbinomial}{}}} -\subsection{Method \code{nbinomial()}}{ +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-monty_rng-negative_binomial_prob}{}}} +\subsection{Method \code{negative_binomial_prob()}}{ Generate \code{n} numbers from a negative binomial distribution \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$nbinomial(n, size, prob, n_threads = 1L)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{monty_rng$negative_binomial_prob(n, size, prob, n_threads = 1L)}\if{html}{\out{
}} } \subsection{Arguments}{ @@ -360,6 +364,31 @@ Generate \code{n} numbers from a negative binomial distribution \item{\code{prob}}{The probability of success on each trial (between 0 and 1, length 1 or n)} +\item{\code{n_threads}}{Number of threads to use; see Details} +} +\if{html}{\out{}} +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-monty_rng-negative_binomial_mu}{}}} +\subsection{Method \code{negative_binomial_mu()}}{ +Generate \code{n} numbers from a negative binomial distribution +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{monty_rng$negative_binomial_mu(n, size, mu, n_threads = 1L)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{n}}{Number of samples to draw (per stream)} + +\item{\code{size}}{The target number of successful trials +(zero or more, length 1 or n)} + +\item{\code{mu}}{The mean +(zero or more, length 1 or n)} + \item{\code{n_threads}}{Number of threads to use; see Details} } \if{html}{\out{
}} diff --git a/src/cpp11.cpp b/src/cpp11.cpp index 7628b53a..b5943761 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -78,10 +78,17 @@ extern "C" SEXP _monty_monty_rng_binomial(SEXP ptr, SEXP n, SEXP r_size, SEXP r_ END_CPP11 } // random.cpp -cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads, bool is_float); -extern "C" SEXP _monty_monty_rng_nbinomial(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP n_threads, SEXP is_float) { +cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads, bool is_float); +extern "C" SEXP _monty_monty_rng_negative_binomial_prob(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP n_threads, SEXP is_float) { BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_nbinomial(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(is_float))); + return cpp11::as_sexp(monty_rng_negative_binomial_prob(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(is_float))); + END_CPP11 +} +// random.cpp +cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_mu, int n_threads, bool is_float); +extern "C" SEXP _monty_monty_rng_negative_binomial_mu(SEXP ptr, SEXP n, SEXP r_size, SEXP r_mu, SEXP n_threads, SEXP is_float) { + BEGIN_CPP11 + return cpp11::as_sexp(monty_rng_negative_binomial_mu(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_mu), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(is_float))); END_CPP11 } // random.cpp @@ -140,6 +147,27 @@ extern "C" SEXP _monty_monty_rng_state(SEXP ptr, SEXP is_float) { return cpp11::as_sexp(monty_rng_state(cpp11::as_cpp>(ptr), cpp11::as_cpp>(is_float))); END_CPP11 } +// random2.cpp +cpp11::doubles cpp_monty_random_real(SEXP ptr); +extern "C" SEXP _monty_cpp_monty_random_real(SEXP ptr) { + BEGIN_CPP11 + return cpp11::as_sexp(cpp_monty_random_real(cpp11::as_cpp>(ptr))); + END_CPP11 +} +// random2.cpp +cpp11::doubles cpp_monty_random_binomial(cpp11::doubles r_size, cpp11::doubles r_prob, SEXP ptr); +extern "C" SEXP _monty_cpp_monty_random_binomial(SEXP r_size, SEXP r_prob, SEXP ptr) { + BEGIN_CPP11 + return cpp11::as_sexp(cpp_monty_random_binomial(cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(ptr))); + END_CPP11 +} +// random2.cpp +cpp11::doubles cpp_monty_random_exponential_rate(cpp11::doubles r_rate, SEXP ptr); +extern "C" SEXP _monty_cpp_monty_random_exponential_rate(SEXP r_rate, SEXP ptr) { + BEGIN_CPP11 + return cpp11::as_sexp(cpp_monty_random_exponential_rate(cpp11::as_cpp>(r_rate), cpp11::as_cpp>(ptr))); + END_CPP11 +} // rng_pointer.cpp cpp11::sexp monty_rng_pointer_init(int n_streams, cpp11::sexp seed, int long_jump, std::string algorithm); extern "C" SEXP _monty_monty_rng_pointer_init(SEXP n_streams, SEXP seed, SEXP long_jump, SEXP algorithm) { @@ -172,29 +200,33 @@ extern "C" SEXP _monty_test_xoshiro_run(SEXP obj) { extern "C" { static const R_CallMethodDef CallEntries[] = { - {"_monty_monty_rng_alloc", (DL_FUNC) &_monty_monty_rng_alloc, 4}, - {"_monty_monty_rng_beta", (DL_FUNC) &_monty_monty_rng_beta, 6}, - {"_monty_monty_rng_binomial", (DL_FUNC) &_monty_monty_rng_binomial, 6}, - {"_monty_monty_rng_cauchy", (DL_FUNC) &_monty_monty_rng_cauchy, 6}, - {"_monty_monty_rng_exponential_mean", (DL_FUNC) &_monty_monty_rng_exponential_mean, 5}, - {"_monty_monty_rng_exponential_rate", (DL_FUNC) &_monty_monty_rng_exponential_rate, 5}, - {"_monty_monty_rng_gamma_rate", (DL_FUNC) &_monty_monty_rng_gamma_rate, 6}, - {"_monty_monty_rng_gamma_scale", (DL_FUNC) &_monty_monty_rng_gamma_scale, 6}, - {"_monty_monty_rng_hypergeometric", (DL_FUNC) &_monty_monty_rng_hypergeometric, 7}, - {"_monty_monty_rng_jump", (DL_FUNC) &_monty_monty_rng_jump, 2}, - {"_monty_monty_rng_long_jump", (DL_FUNC) &_monty_monty_rng_long_jump, 2}, - {"_monty_monty_rng_multinomial", (DL_FUNC) &_monty_monty_rng_multinomial, 6}, - {"_monty_monty_rng_nbinomial", (DL_FUNC) &_monty_monty_rng_nbinomial, 6}, - {"_monty_monty_rng_normal", (DL_FUNC) &_monty_monty_rng_normal, 7}, - {"_monty_monty_rng_pointer_init", (DL_FUNC) &_monty_monty_rng_pointer_init, 4}, - {"_monty_monty_rng_pointer_sync", (DL_FUNC) &_monty_monty_rng_pointer_sync, 2}, - {"_monty_monty_rng_poisson", (DL_FUNC) &_monty_monty_rng_poisson, 5}, - {"_monty_monty_rng_random_normal", (DL_FUNC) &_monty_monty_rng_random_normal, 5}, - {"_monty_monty_rng_random_real", (DL_FUNC) &_monty_monty_rng_random_real, 4}, - {"_monty_monty_rng_state", (DL_FUNC) &_monty_monty_rng_state, 2}, - {"_monty_monty_rng_uniform", (DL_FUNC) &_monty_monty_rng_uniform, 6}, - {"_monty_test_rng_pointer_get", (DL_FUNC) &_monty_test_rng_pointer_get, 2}, - {"_monty_test_xoshiro_run", (DL_FUNC) &_monty_test_xoshiro_run, 1}, + {"_monty_cpp_monty_random_binomial", (DL_FUNC) &_monty_cpp_monty_random_binomial, 3}, + {"_monty_cpp_monty_random_exponential_rate", (DL_FUNC) &_monty_cpp_monty_random_exponential_rate, 2}, + {"_monty_cpp_monty_random_real", (DL_FUNC) &_monty_cpp_monty_random_real, 1}, + {"_monty_monty_rng_alloc", (DL_FUNC) &_monty_monty_rng_alloc, 4}, + {"_monty_monty_rng_beta", (DL_FUNC) &_monty_monty_rng_beta, 6}, + {"_monty_monty_rng_binomial", (DL_FUNC) &_monty_monty_rng_binomial, 6}, + {"_monty_monty_rng_cauchy", (DL_FUNC) &_monty_monty_rng_cauchy, 6}, + {"_monty_monty_rng_exponential_mean", (DL_FUNC) &_monty_monty_rng_exponential_mean, 5}, + {"_monty_monty_rng_exponential_rate", (DL_FUNC) &_monty_monty_rng_exponential_rate, 5}, + {"_monty_monty_rng_gamma_rate", (DL_FUNC) &_monty_monty_rng_gamma_rate, 6}, + {"_monty_monty_rng_gamma_scale", (DL_FUNC) &_monty_monty_rng_gamma_scale, 6}, + {"_monty_monty_rng_hypergeometric", (DL_FUNC) &_monty_monty_rng_hypergeometric, 7}, + {"_monty_monty_rng_jump", (DL_FUNC) &_monty_monty_rng_jump, 2}, + {"_monty_monty_rng_long_jump", (DL_FUNC) &_monty_monty_rng_long_jump, 2}, + {"_monty_monty_rng_multinomial", (DL_FUNC) &_monty_monty_rng_multinomial, 6}, + {"_monty_monty_rng_negative_binomial_mu", (DL_FUNC) &_monty_monty_rng_negative_binomial_mu, 6}, + {"_monty_monty_rng_negative_binomial_prob", (DL_FUNC) &_monty_monty_rng_negative_binomial_prob, 6}, + {"_monty_monty_rng_normal", (DL_FUNC) &_monty_monty_rng_normal, 7}, + {"_monty_monty_rng_pointer_init", (DL_FUNC) &_monty_monty_rng_pointer_init, 4}, + {"_monty_monty_rng_pointer_sync", (DL_FUNC) &_monty_monty_rng_pointer_sync, 2}, + {"_monty_monty_rng_poisson", (DL_FUNC) &_monty_monty_rng_poisson, 5}, + {"_monty_monty_rng_random_normal", (DL_FUNC) &_monty_monty_rng_random_normal, 5}, + {"_monty_monty_rng_random_real", (DL_FUNC) &_monty_monty_rng_random_real, 4}, + {"_monty_monty_rng_state", (DL_FUNC) &_monty_monty_rng_state, 2}, + {"_monty_monty_rng_uniform", (DL_FUNC) &_monty_monty_rng_uniform, 6}, + {"_monty_test_rng_pointer_get", (DL_FUNC) &_monty_test_rng_pointer_get, 2}, + {"_monty_test_xoshiro_run", (DL_FUNC) &_monty_test_xoshiro_run, 1}, {NULL, NULL, 0} }; } diff --git a/src/random.cpp b/src/random.cpp index cc16b0e4..88e95cf8 100644 --- a/src/random.cpp +++ b/src/random.cpp @@ -336,7 +336,7 @@ cpp11::sexp monty_rng_binomial(SEXP ptr, int n, template -cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, +cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads) { T *rng = cpp11::as_cpp>(ptr).get(); @@ -363,7 +363,7 @@ cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, for (size_t j = 0; j < (size_t)n; ++j) { auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; auto prob_ij = prob_vary.draw ? prob_i[j] : prob_i[0]; - y_i[j] = monty::random::nbinomial(state, size_ij, prob_ij); + y_i[j] = monty::random::negative_binomial_prob(state, size_ij, prob_ij); } } catch (std::exception const& e) { errors.capture(e, i); @@ -375,6 +375,48 @@ cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, return sexp_matrix(ret, n, n_streams); } + +template +cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, + cpp11::doubles r_size, cpp11::doubles r_mu, + int n_threads) { + T *rng = cpp11::as_cpp>(ptr).get(); + const int n_streams = rng->size(); + cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); + double * y = REAL(ret); + + const double * size = REAL(r_size); + const double * mu = REAL(r_mu); + auto size_vary = check_input_type(r_size, n, n_streams, "size"); + auto mu_vary = check_input_type(r_mu, n, n_streams, "mu"); + + monty::utils::openmp_errors errors(n_streams); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) num_threads(n_threads) +#endif + for (int i = 0; i < n_streams; ++i) { + try { + auto &state = rng->state(i); + auto y_i = y + n * i; + auto size_i = size_vary.generator ? size + size_vary.offset * i : size; + auto mu_i = mu_vary.generator ? mu + mu_vary.offset * i : mu; + for (size_t j = 0; j < (size_t)n; ++j) { + auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; + auto mu_ij = mu_vary.draw ? mu_i[j] : mu_i[0]; + y_i[j] = monty::random::negative_binomial_mu(state, size_ij, mu_ij); + } + } catch (std::exception const& e) { + errors.capture(e, i); + } + } + + errors.report("generators", 4, true); + + return sexp_matrix(ret, n, n_streams); +} + + template cpp11::sexp monty_rng_poisson(SEXP ptr, int n, cpp11::doubles r_lambda, int n_threads) { @@ -812,12 +854,21 @@ cpp11::sexp monty_rng_binomial(SEXP ptr, int n, } [[cpp11::register]] -cpp11::sexp monty_rng_nbinomial(SEXP ptr, int n, +cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads, bool is_float) { return is_float ? - monty_rng_nbinomial(ptr, n, r_size, r_prob, n_threads) : - monty_rng_nbinomial(ptr, n, r_size, r_prob, n_threads); + monty_rng_negative_binomial_prob(ptr, n, r_size, r_prob, n_threads) : + monty_rng_negative_binomial_prob(ptr, n, r_size, r_prob, n_threads); +} + +[[cpp11::register]] +cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, + cpp11::doubles r_size, cpp11::doubles r_mu, + int n_threads, bool is_float) { + return is_float ? + monty_rng_negative_binomial_mu(ptr, n, r_size, r_mu, n_threads) : + monty_rng_negative_binomial_mu(ptr, n, r_size, r_mu, n_threads); } [[cpp11::register]] diff --git a/src/random2.cpp b/src/random2.cpp new file mode 100644 index 00000000..5786f0e4 --- /dev/null +++ b/src/random2.cpp @@ -0,0 +1,73 @@ +#include +#include +#include +#include + +#include +#include +#include + +using default_rng64 = monty::random::prng>; + +template +void check_length(T x, size_t len, const char * name) { + const size_t len_given = Rf_length(x); + if (len_given != len) { + cpp11::stop("Expected '%s' to have length %d, not %d", + name, len, len_given); + } +} + + +[[cpp11::register]] +cpp11::doubles cpp_monty_random_real(SEXP ptr) { + default_rng64 *rng = cpp11::as_cpp>(ptr).get(); + const int n = rng->size(); + cpp11::writable::doubles r_y = cpp11::writable::doubles(n); + double *y = REAL(r_y); + for (int i = 0; i < n; ++i) { + y[i] = monty::random::random_real(rng->state(i)); + } + return r_y; +} + + +[[cpp11::register]] +cpp11::doubles cpp_monty_random_binomial(cpp11::doubles r_size, + cpp11::doubles r_prob, + SEXP ptr) { + default_rng64 *rng = cpp11::as_cpp>(ptr).get(); + const int n = rng->size(); + check_length(r_size, n, "size"); + check_length(r_prob, n, "size"); + const double * size = REAL(r_size); + const double * prob = REAL(r_prob); + + cpp11::writable::doubles r_y = cpp11::writable::doubles(n); + double *y = REAL(r_y); + for (int i = 0; i < n; ++i) { + auto &state = rng->state(i); + y[i] = monty::random::binomial(state, size[i], prob[i]); + } + + return r_y; +} + + +[[cpp11::register]] +cpp11::doubles cpp_monty_random_exponential_rate(cpp11::doubles r_rate, + SEXP ptr) { + default_rng64 *rng = cpp11::as_cpp>(ptr).get(); + const int n = rng->size(); + check_length(r_rate, n, "size"); + const double * rate = REAL(r_rate); + + cpp11::writable::doubles r_y = cpp11::writable::doubles(n); + double *y = REAL(r_y); + for (int i = 0; i < n; ++i) { + auto &state = rng->state(i); + y[i] = monty::random::exponential_rate(state, rate[i]); + } + + return r_y; +} diff --git a/tests/testthat/helper-monty.R b/tests/testthat/helper-monty.R index ae4982b6..5308e6da 100644 --- a/tests/testthat/helper-monty.R +++ b/tests/testthat/helper-monty.R @@ -81,154 +81,6 @@ ex_simple_nested_with_base <- function(n_groups) { } -ex_simple_nested_with_base_stochastic <- function(n_groups) { - e <- new.env(parent = topenv()) - e$sigma <- 5 - e$mu <- rnorm(n_groups, 0, e$sigma) - e$n_groups <- n_groups - with( - e, - monty_model(list( - parameters = c("sigma", paste0("mu_", seq_len(n_groups))), - direct_sample = function(rng) { - sigma <- rng$uniform(1, 0, 10) - alpha <- rng$uniform(n_groups, 0, 1) - c(sigma, rng$normal(n_groups, 0, alpha * sigma)) - }, - density = function(x, by_group = FALSE) { - density1 <- function(y) { - sigma <- y[[1]] - if (sigma <= 0) { - rep(-Inf, length(y) - 1) - } else { - alpha <- runif(1, 0, 1) - dnorm(y[-1], 0, alpha * y[[1]], log = TRUE) - } - } - if (is.matrix(x)) { - z <- vapply(seq_len(ncol(x)), function(i) density1(x[, i]), - numeric(nrow(x) - 1)) - value <- colSums(z) + dunif(x[1, ], 0, 10, log = TRUE) - if (by_group) structure(value, "by_group" = z) else value - } else { - z <- density1(x) - value <- sum(z) + dunif(x[[1]], 0, 10, log = TRUE) - if (by_group) structure(value, "by_group" = z) else value - } - }, - parameter_groups = c(0, seq_len(n_groups)), - mu = mu), - monty_model_properties(allow_multiple_parameters = TRUE, - is_stochastic = TRUE))) -} - - -ex_dust_sir <- function(n_particles = 100, n_threads = 1, - deterministic = FALSE, save_trajectories = FALSE) { - testthat::skip_if_not_installed("dust") - sir <- dust::dust_example("sir") - - np <- 10 - end <- 150 * 4 - times <- seq(0, end, by = 4) - ans <- sir$new(list(), 0, np, seed = 1L)$simulate(times) - dat <- data.frame(time = times[-1], incidence = ans[5, 1, -1]) - - ## TODO: an upshot here is that our dust models are always going to - ## need to be initialisable; we might need to sample from the - ## statistical parameters, or set things up to allow two-phases of - ## initialsation (which is I think where we are heading, so that's - ## fine). - model <- sir$new(list(), 0, n_particles, seed = 1L, n_threads = n_threads, - deterministic = deterministic) - model$set_data(dust::dust_data(dat)) - model$set_index(c(2, 4)) - - prior_beta_shape <- 1 - prior_beta_rate <- 1 / 0.5 - prior_gamma_shape <- 1 - prior_gamma_rate <- 1 / 0.5 - - trajectories <- NULL - - density <- function(x) { - beta <- x[[1]] - gamma <- x[[2]] - prior <- dgamma(beta, prior_beta_shape, prior_beta_rate, log = TRUE) + - dgamma(gamma, prior_gamma_shape, prior_gamma_rate, log = TRUE) - if (is.finite(prior)) { - model$update_state( - pars = list(beta = x[[1]], gamma = x[[2]]), - time = 0, - set_initial_state = TRUE) - res <- model$filter(save_trajectories = save_trajectories) - if (save_trajectories) { - trajectories <<- res$trajectories - } - ll <- res$log_likelihood - } else { - ll <- -Inf - } - ll + prior - } - - direct_sample <- function(rng) { - c(rng$gamma_scale(1, prior_beta_shape, 1 / prior_beta_rate), - rng$gamma_scale(1, prior_gamma_shape, 1 / prior_gamma_rate)) - } - - set_rng_state <- function(rng_state) { - n_streams <- n_particles + 1 - if (length(rng_state) != 32 * n_streams) { - ## Expand the state by short jumps; we'll make this nicer once - ## we refactor the RNG interface and dust. - rng_state <- monty_rng$new(rng_state, n_streams)$state() - } - model$set_rng_state(rng_state) - } - - get_rng_state <- function() { - model$rng_state() - } - - if (save_trajectories) { - observer <- monty_observer( - function() { - ## TODO: It's not really clear to me (Rich) that we want the - ## rng coming in here. In dust2 we'll correctly use the - ## *filter* rng. however, even there we end up with the act - ## of observation changing the behaviour of the sampler and I - ## am not sure that is really desirable. We could probably - ## improve that in dust, but that would require that we do not - ## pass the rng here too. So for now we take the first - ## particle. - ## i <- floor(rng$random_real(1) * model$model$n_particles()) + 1L - i <- 1L - if (save_trajectories) { - traj <- trajectories[, i, , drop = FALSE] - dim(traj) <- dim(traj)[-2] - } else { - traj <- NULL - } - list(trajectories = traj, state = model$state()[, i]) - }) - } else { - observer <- NULL - } - - monty_model( - list(model = model, - density = density, - direct_sample = direct_sample, - parameters = c("beta", "gamma"), - domain = cbind(c(0, 0), c(Inf, Inf)), - observer = observer, - set_rng_state = set_rng_state, - get_rng_state = get_rng_state), - monty_model_properties(is_stochastic = !deterministic)) -} - - random_array <- function(dim, named = FALSE) { if (named) { dn <- lapply(seq_along(dim), function(i) { @@ -242,91 +94,21 @@ random_array <- function(dim, named = FALSE) { } -ex_dust_sir_likelihood <- function(n_particles = 100, n_threads = 1, - deterministic = FALSE, - save_trajectories = FALSE) { - testthat::skip_if_not_installed("dust") - sir <- dust::dust_example("sir") - - np <- 10 - end <- 150 * 4 - times <- seq(0, end, by = 4) - ans <- sir$new(list(), 0, np, seed = 1L)$simulate(times) - dat <- data.frame(time = times[-1], incidence = ans[5, 1, -1]) - - ## TODO: an upshot here is that our dust models are always going to - ## need to be initialisable; we might need to sample from the - ## statistical parameters, or set things up to allow two-phases of - ## initialsation (which is I think where we are heading, so that's - ## fine). - model <- sir$new(list(), 0, n_particles, seed = 1L, n_threads = n_threads, - deterministic = deterministic) - model$set_data(dust::dust_data(dat)) - model$set_index(c(2, 4)) - - trajectories <- NULL - - density <- function(x) { - beta <- x[[1]] - gamma <- x[[2]] - model$update_state( - pars = list(beta = x[[1]], gamma = x[[2]]), - time = 0, - set_initial_state = TRUE) - res <- model$filter(save_trajectories = save_trajectories) - if (save_trajectories) { - trajectories <<- res$trajectories - } - res$log_likelihood - } - - set_rng_state <- function(rng_state) { - n_streams <- n_particles + 1 - if (length(rng_state) != 32 * n_streams) { - ## Expand the state by short jumps; we'll make this nicer once - ## we refactor the RNG interface and dust. - rng_state <- monty_rng$new(rng_state, n_streams)$state() - } - model$set_rng_state(rng_state) - } - - get_rng_state <- function() { - model$rng_state() - } +ex_sir_filter_likelihood <- function(n_particles = 100, + deterministic = FALSE, + save_trajectories = FALSE) { + data <- data.frame(time = c( 4, 8, 12, 16, 20, 24, 28, 32, 36), + incidence = c( 1, 0, 3, 5, 2, 4, 3, 7, 2)) + sir_filter_monty(data, n_particles, deterministic, save_trajectories) +} - if (save_trajectories) { - observer <- monty_observer( - function() { - ## TODO: It's not really clear to me (Rich) that we want the - ## rng coming in here. In dust2 we'll correctly use the - ## *filter* rng. however, even there we end up with the act - ## of observation changing the behaviour of the sampler and I - ## am not sure that is really desirable. We could probably - ## improve that in dust, but that would require that we do not - ## pass the rng here too. So for now we take the first - ## particle. - ## i <- floor(rng$random_real(1) * model$model$n_particles()) + 1L - i <- 1L - if (save_trajectories) { - traj <- trajectories[, i, , drop = FALSE] - dim(traj) <- dim(traj)[-2] - } else { - traj <- NULL - } - list(trajectories = traj, state = model$state()[, i]) - }) - } else { - observer <- NULL - } - monty_model( - list(model = model, - density = density, - parameters = c("beta", "gamma"), - observer = observer, - set_rng_state = set_rng_state, - get_rng_state = get_rng_state), - monty_model_properties(is_stochastic = !deterministic)) +ex_sir_filter_posterior <- function(...) { + prior <- monty_dsl({ + beta ~ Gamma(shape = 1, rate = 1 / 0.5) + gamma ~ Gamma(shape = 1, rate = 1 / 0.5) + }) + ex_sir_filter_likelihood(...) + prior } diff --git a/tests/testthat/helper-sir-filter.R b/tests/testthat/helper-sir-filter.R new file mode 100644 index 00000000..c5cce7a5 --- /dev/null +++ b/tests/testthat/helper-sir-filter.R @@ -0,0 +1,123 @@ +sir_filter_monty <- function(data, n_particles, deterministic = FALSE, + save_trajectories = FALSE, seed = NULL) { + parameters <- c("beta", "gamma") + env <- new.env() + base <- list(N = 1000, I0 = 10, beta = 0.2, gamma = 0.1, exp_noise = 1e6) + + get_rng_state <- function() { + c(monty_rng_state(env$rng$filter$ptr, FALSE), + monty_rng_state(env$rng$system$ptr, FALSE)) + } + + set_rng_state <- function(rng_state) { + n_streams <- n_particles + 1 + r <- matrix( + monty_rng$new(n_streams = n_streams, seed = rng_state)$state(), + ncol = n_streams) + env$rng <- list( + filter = monty_random_alloc(1, r[, 1], deterministic), + system = monty_random_alloc(n_particles, c(r[, -1]), deterministic)) + } + + set_rng_state(seed) + + density <- function(x) { + pars <- base + pars[parameters] <- x + res <- sir_filter(pars, data, n_particles, env$rng, save_trajectories) + if (save_trajectories) { + env$trajectories <- res$trajectories + } + res$log_likelihood + } + + if (save_trajectories) { + observer <- monty_observer( + function() { + i <- 1L + trajectories <- env$trajectories[c(2, 4), i, , drop = FALSE] + dim(trajectories) <- dim(trajectories)[-2] + list(trajectories = trajectories) + }) + } else { + observer <- NULL + } + + monty_model( + list(density = density, + parameters = parameters, + observer = observer, + set_rng_state = set_rng_state, + get_rng_state = get_rng_state), + monty_model_properties(is_stochastic = !deterministic)) +} + + +sir_filter <- function(pars, data, n_particles, rng, + save_trajectories = FALSE) { + y <- list(S = pars$N - pars$I0, I = pars$I0, R = 0, cases = 0) + packer <- monty_packer(names(y)) + state <- matrix(packer$pack(y), 4, n_particles) + time <- 0 + dt <- 1 + ll <- 0 + exp_noise <- rep_len(pars$exp_noise, n_particles) + i_cases <- packer$index()$cases + if (save_trajectories) { + trajectories <- array(NA_real_, c(length(y), n_particles, nrow(data))) + } else { + trajectories <- NULL + } + + for (i in seq_len(nrow(data))) { + from <- time + to <- data$time[[i]] + state <- sir_run(from, to, dt, state, packer, pars, rng$system) + noise <- monty_random_exponential_rate(exp_noise, rng$system) + lambda <- state[i_cases, ] + noise + tmp <- dpois(data$incidence[[i]], lambda, log = TRUE) + w <- exp(tmp - max(tmp)) + ll <- ll + log(mean(w)) + max(tmp) + u <- monty_random_real(rng$filter) + k <- dust_resample_weight(w, u) + + if (save_trajectories) { + trajectories[, , i] <- state + trajectories <- trajectories[, k, , drop = FALSE] # slow but correct... + } + + state <- state[, k, drop = FALSE] + time <- to + } + + list(log_likelihood = ll, trajectories = trajectories) +} + + +sir_run <- function(from, to, dt, state, packer, pars, rng) { + state <- packer$unpack(state) + for (time in seq(from, to, by = dt)[-1]) { + state <- sir_step(time, dt, state, pars, rng) + } + packer$pack(state) +} + + +sir_step <- function(time, dt, state, pars, rng) { + p_SI <- 1 - exp(-pars$beta * state$I / pars$N * dt) + p_IR <- rep(1 - exp(-pars$gamma * dt), length(state$I)) + n_SI <- monty_random_binomial(state$S, p_SI, rng) + n_IR <- monty_random_binomial(state$I, p_IR, rng) + cases <- if (time %% 1 == 0) 0 else state$cases + list(S = state$S - n_SI, + I = state$I + n_SI - n_IR, + R = state$R + n_IR, + cases = cases + n_SI) +} + + +dust_resample_weight <- function(w, u) { + n <- length(w) + uu <- u / n + seq(0, by = 1 / n, length.out = n) + findInterval(uu, cumsum(w / sum(w))) + 1L +} diff --git a/tests/testthat/test-combine.R b/tests/testthat/test-combine.R index 525415c0..21fa709e 100644 --- a/tests/testthat/test-combine.R +++ b/tests/testthat/test-combine.R @@ -186,7 +186,7 @@ test_that("can combine gradients where parameters do not agree", { test_that("can combine a stochastic and deterministic model", { - ll <- ex_dust_sir_likelihood() + ll <- ex_sir_filter_likelihood() prior <- monty_dsl({ beta ~ Gamma(shape = 1, rate = 1 / 0.5) gamma ~ Gamma(shape = 1, rate = 1 / 0.5) @@ -199,7 +199,7 @@ test_that("can combine a stochastic and deterministic model", { test_that("can't disable stochastic model on combination", { - ll <- ex_dust_sir_likelihood() + ll <- ex_sir_filter_likelihood() prior <- monty_dsl({ beta ~ Gamma(shape = 1, rate = 1 / 0.5) gamma ~ Gamma(shape = 1, rate = 1 / 0.5) @@ -212,7 +212,7 @@ test_that("can't disable stochastic model on combination", { test_that("can't create model out of two stochastic halves", { - ll <- ex_dust_sir_likelihood() + ll <- ex_sir_filter_likelihood() expect_error( ll + ll, "Can't combine two stochastic models") diff --git a/tests/testthat/test-packer.R b/tests/testthat/test-packer.R index 5025b1d6..ec688ac2 100644 --- a/tests/testthat/test-packer.R +++ b/tests/testthat/test-packer.R @@ -6,7 +6,7 @@ test_that("can't create empty packer", { test_that("trivial packer", { xp <- monty_packer("a") - expect_equal(xp$parameters, "a") + expect_equal(xp$names(), "a") expect_equal(xp$unpack(1), list(a = 1)) expect_equal(xp$unpack(c(a = 1)), list(a = 1)) expect_error(xp$unpack(c(b = 1)), @@ -23,7 +23,7 @@ test_that("trivial packer", { test_that("multiple scalar unpacking", { xp <- monty_packer(c("a", "b", "c")) - expect_equal(xp$parameters, c("a", "b", "c")) + expect_equal(xp$names(), c("a", "b", "c")) expect_equal(xp$unpack(1:3), list(a = 1, b = 2, c = 3)) expect_equal(xp$pack(list(a = 1, b = 2, c = 3)), 1:3) expect_equal(xp$index(), list(a = 1, b = 2, c = 3)) @@ -32,7 +32,7 @@ test_that("multiple scalar unpacking", { test_that("can bind data into an unpacked list", { xp <- monty_packer(c("a", "b"), fixed = list(x = 1:5, y = 10)) - expect_equal(xp$parameters, c("a", "b")) + expect_equal(xp$names(), c("a", "b")) expect_equal(xp$unpack(1:2), list(a = 1, b = 2, x = 1:5, y = 10)) expect_equal(xp$pack(list(a = 1, b = 2, x = 1:5, y = 10)), 1:2) expect_equal(xp$pack(list(a = 1, b = 2)), 1:2) @@ -41,14 +41,14 @@ test_that("can bind data into an unpacked list", { test_that("can unpack arrays", { xp <- monty_packer("a", list(b = 3)) - expect_equal(xp$parameters, c("a", "b[1]", "b[2]", "b[3]")) + expect_equal(xp$names(), c("a", "b[1]", "b[2]", "b[3]")) expect_equal(xp$unpack(1:4), list(a = 1, b = 2:4)) }) test_that("can use integer vectors for array inputs", { xp <- monty_packer("a", c(b = 3, c = 4)) - expect_equal(xp$parameters, + expect_equal(xp$names(), c("a", sprintf("b[%d]", 1:3), sprintf("c[%d]", 1:4))) expect_equal(xp$unpack(1:8), list(a = 1, b = 2:4, c = 5:8)) expect_equal(xp$index(), list(a = 1, b = 2:4, c = 5:8)) @@ -65,7 +65,7 @@ test_that("can create packers with higher-level dimensionsality", { xp <- monty_packer( array = list(a = 1, b = 2, c = 2:3, d = 2:4)) expect_equal( - xp$parameters, + xp$names(), c("a[1]", "b[1]", "b[2]", sprintf("c[%d,%d]", 1:2, rep(1:3, each = 2)), sprintf("d[%d,%d,%d]", 1:2, rep(1:3, each = 2), rep(1:4, each = 6)))) @@ -112,7 +112,7 @@ test_that("validate array inputs", { test_that("can pass empty array elements as scalars", { p <- monty_packer(array = list(a = integer(), b = 2)) - expect_equal(p$parameters, c("a", "b[1]", "b[2]")) + expect_equal(p$names(), c("a", "b[1]", "b[2]")) expect_equal(p$unpack(1:3), list(a = 1, b = 2:3)) expect_equal(p$pack(list(a = 1, b = 2:3)), 1:3) }) @@ -120,7 +120,7 @@ test_that("can pass empty array elements as scalars", { test_that("can pass empty array elements as scalars in odd order", { p <- monty_packer(array = list(a = 2, b = NULL)) - expect_equal(p$parameters, c("a[1]", "a[2]", "b")) + expect_equal(p$names(), c("a[1]", "a[2]", "b")) expect_equal(p$unpack(1:3), list(a = 1:2, b = 3)) expect_equal(p$pack(list(a = 1:2, b = 3)), 1:3) }) @@ -138,7 +138,7 @@ test_that("can post-process parameters", { list(d = x$a + x$b + x$c) } xp <- monty_packer(c("a", "b", "c"), process = p) - expect_equal(xp$parameters, c("a", "b", "c")) + expect_equal(xp$names(), c("a", "b", "c")) expect_equal(xp$unpack(1:3), list(a = 1, b = 2, c = 3, d = 6)) expect_equal(xp$pack(xp$unpack(1:3)), 1:3) expect_equal(xp$pack(list(a = 1, b = 2, c = 3)), 1:3) @@ -158,7 +158,7 @@ test_that("require that process is well-behaved", { } xp <- monty_packer(c("a", "b", "c"), process = p) expect_error(xp$unpack(1:3), - "'process()' is trying to overwrite entries in parameters", + "'process()' is trying to overwrite entries in your list", fixed = TRUE) }) @@ -388,7 +388,7 @@ test_that("can subset a packer of scalars", { p <- monty_packer(c("a", "b", "c", "d")) res <- p$subset(c("b", "c")) expect_equal(res$index, 2:3) - expect_equal(res$packer$parameters, c("b", "c")) + expect_equal(res$packer$names(), c("b", "c")) expect_equal(res$packer$unpack(1:2), list(b = 1, c = 2)) }) @@ -399,7 +399,7 @@ test_that("can subset a packer of arrays", { expect_equal(res$index, c(1, 4:12)) cmp <- monty_packer(array = list(a = integer(), c = c(3, 3))) - expect_equal(res$packer$parameters, cmp$parameters) + expect_equal(res$packer$names(), cmp$names()) expect_equal(res$packer$index(), cmp$index()) }) diff --git a/tests/testthat/test-rng.R b/tests/testthat/test-rng.R index 41de599b..2b022aa4 100644 --- a/tests/testthat/test-rng.R +++ b/tests/testthat/test-rng.R @@ -1357,13 +1357,24 @@ test_that("can generate negative binomial numbers", { m <- 1000000 n <- 958 p <- 0.004145 - yf <- monty_rng$new(1)$nbinomial(m, n, p) + yf <- monty_rng$new(1)$negative_binomial_prob(m, n, p) expect_equal(mean(yf), (1 - p) * n / p, tolerance = 1e-3) expect_equal(var(yf), ((1 - p) * n) / p^2, tolerance = 1e-2) }) +test_that("negative_binomial_mu follows from negative_binomial_prob", { + rng1 <- monty_rng$new(seed = 1L) + rng2 <- monty_rng$new(seed = 1L) + prob <- 0.3 + size <- 20 + expect_identical( + rng1$negative_binomial_mu(100, size, size * (1 - prob) / prob), + rng2$negative_binomial_prob(100, size, prob)) +}) + + test_that("deterministic negative binomial returns mean", { m <- 100 p <- as.numeric(sample(10, m, replace = TRUE)) / 10 @@ -1372,22 +1383,23 @@ test_that("deterministic negative binomial returns mean", { rng_f <- monty_rng$new(1, real_type = "float", deterministic = TRUE) rng_d <- monty_rng$new(1, real_type = "double", deterministic = TRUE) - expect_equal(rng_f$nbinomial(m, n, p), (1 - p) * n / p, tolerance = 1e-6) - expect_equal(rng_d$nbinomial(m, n, p), (1 - p) * n / p) + expect_equal(rng_f$negative_binomial_prob(m, n, p), (1 - p) * n / p, + tolerance = 1e-6) + expect_equal(rng_d$negative_binomial_prob(m, n, p), (1 - p) * n / p) }) test_that("negative binomial prevents bad inputs", { - expect_error(monty_rng$new(1)$nbinomial(1, 10, 0), - "Invalid call to nbinomial with size = 10, prob = 0") - expect_error(monty_rng$new(1)$nbinomial(1, 0, 0.5), - "Invalid call to nbinomial with size = 0, prob = 0.5") - expect_error(monty_rng$new(1)$nbinomial(1, 10, 1.5), - "Invalid call to nbinomial with size = 10, prob = 1.5") - expect_error(monty_rng$new(1)$nbinomial(1, 10, Inf), - "Invalid call to nbinomial with size = 10, prob = inf") - expect_error(monty_rng$new(1)$nbinomial(1, Inf, 0.4), - "Invalid call to nbinomial with size = inf, prob = 0.4") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 10, 0), + "Invalid call to negative_binomial with size = 10, prob = 0") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 0, 0.5), + "Invalid call to negative_binomial with size = 0, prob = 0.5") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 10, 1.5), + "Invalid call to negative_binomial with size = 10, prob = 1.5") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, 10, Inf), + "Invalid call to negative_binomial with size = 10, prob = inf") + expect_error(monty_rng$new(1)$negative_binomial_prob(1, Inf, 0.4), + "Invalid call to negative_binomial with size = inf, prob = 0.4") }) diff --git a/tests/testthat/test-sample.R b/tests/testthat/test-sample.R index cdb381fe..82117b2e 100644 --- a/tests/testthat/test-sample.R +++ b/tests/testthat/test-sample.R @@ -269,7 +269,7 @@ test_that("error if initial conditions do not have finite density", { test_that("can continue a stochastic model identically", { set.seed(1) - model <- ex_dust_sir() + model <- ex_sir_filter_posterior() vcv <- matrix(c(0.0006405, 0.0005628, 0.0005628, 0.0006641), 2, 2) sampler <- monty_sampler_random_walk(vcv = vcv) initial <- c(0.2, 0.1) diff --git a/tests/testthat/test-sampler-adaptive.R b/tests/testthat/test-sampler-adaptive.R index 249348de..3b4bce4c 100644 --- a/tests/testthat/test-sampler-adaptive.R +++ b/tests/testthat/test-sampler-adaptive.R @@ -92,7 +92,7 @@ test_that("can continue adaptive sampler", { test_that("can't use adaptive sampler with stochastic models", { set.seed(1) - m <- ex_dust_sir() + m <- ex_sir_filter_posterior() sampler <- monty_sampler_adaptive(initial_vcv = diag(c(0.01, 0.01))) expect_error( monty_sample(m, sampler, 30, n_chains = 3), diff --git a/tests/testthat/test-sampler-hmc.R b/tests/testthat/test-sampler-hmc.R index b39b13f2..0383b41c 100644 --- a/tests/testthat/test-sampler-hmc.R +++ b/tests/testthat/test-sampler-hmc.R @@ -216,7 +216,7 @@ test_that("can't use hmc with models that lack gradients", { test_that("can't use hmc with stochastic models", { set.seed(1) - m <- ex_dust_sir() + m <- ex_sir_filter_posterior() sampler <- monty_sampler_hmc(epsilon = 0.1, n_integration_steps = 10) expect_error( monty_sample(m, sampler, 30, n_chains = 3), diff --git a/tests/testthat/test-sampler-random-walk.R b/tests/testthat/test-sampler-random-walk.R index 88e19e8e..8229b709 100644 --- a/tests/testthat/test-sampler-random-walk.R +++ b/tests/testthat/test-sampler-random-walk.R @@ -26,7 +26,7 @@ test_that("validate sampler against model on initialisation", { test_that("can draw samples from a random model", { set.seed(1) - m <- ex_dust_sir() + m <- ex_sir_filter_posterior() vcv <- matrix(c(0.0006405, 0.0005628, 0.0005628, 0.0006641), 2, 2) sampler <- monty_sampler_random_walk(vcv = vcv) res <- monty_sample(m, sampler, 20) @@ -36,7 +36,7 @@ test_that("can draw samples from a random model", { test_that("can observe a model", { - m <- ex_dust_sir(save_trajectories = TRUE) + m <- ex_sir_filter_posterior(save_trajectories = TRUE) vcv <- matrix(c(0.0006405, 0.0005628, 0.0005628, 0.0006641), 2, 2) sampler <- monty_sampler_random_walk(vcv = vcv) @@ -46,18 +46,14 @@ test_that("can observe a model", { expect_setequal(names(res), c("pars", "density", "initial", "details", "observations")) expect_equal(names(res$observations), - c("trajectories", "state")) + "trajectories") expect_equal(dim(res$observations$trajectories), - c(2, 151, 20, 3)) # states, time steps, samples, chains - expect_equal(dim(res$observations$state), - c(5, 20, 3)) # states, samples, chains - expect_equal(res$observations$state[c(2, 4), , ], - res$observations$trajectories[, 151, , ]) + c(2, 9, 20, 3)) # states, time steps, samples, chains }) test_that("can continue observed models", { - m <- ex_dust_sir(save_trajectories = TRUE) + m <- ex_sir_filter_posterior(save_trajectories = TRUE) vcv <- matrix(c(0.0006405, 0.0005628, 0.0005628, 0.0006641), 2, 2) sampler <- monty_sampler_random_walk(vcv = vcv) @@ -332,7 +328,7 @@ test_that("can print a sampler object", { test_that("Can rerun a stochastic model", { set.seed(1) - m <- ex_dust_sir() + m <- ex_sir_filter_posterior() vcv <- matrix(c(0.0006405, 0.0005628, 0.0005628, 0.0006641), 2, 2) sampler1 <- monty_sampler_random_walk(vcv = vcv, rerun_every = 2) sampler2 <- monty_sampler_random_walk(vcv = vcv, rerun_every = 2,