From 14b075b044cead567e19f4857e48621c396bc527 Mon Sep 17 00:00:00 2001 From: edknock Date: Thu, 5 Sep 2024 16:44:47 +0100 Subject: [PATCH 01/14] implement rerunfor monty_sampler_random_walk --- R/sampler-random-walk.R | 51 +++++++++++++++++++++++++++++++- man/monty_sampler_random_walk.Rd | 23 +++++++++++++- 2 files changed, 72 insertions(+), 2 deletions(-) diff --git a/R/sampler-random-walk.R b/R/sampler-random-walk.R index abb909df..09310900 100644 --- a/R/sampler-random-walk.R +++ b/R/sampler-random-walk.R @@ -22,16 +22,38 @@ ##' ##' The initial point selected will lie within the domain, as this is ##' enforced by [monty_sample]. +##' +##' @param 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 +##' (`Inf`) will never rerun this point, but if you set to 100, then +##' every 100 steps we run the model on both the proposed *and* previously +##' accepted parameters before doing the comparison. This may help "unstick" +##' chains, at the cost of some bias in the results. +##' +##' @param rerun_random Logical, controlling the behaviour of +##' rerunning (when `rerun_every` is finite). The default value of +##' `FALSE` will rerun the model at fixed intervals of iterations +##' (given by `rerun_every`). If `TRUE`, then we stochastically +##' rerun each step with probability of `1 / rerun_every`. This gives +##' the same expected number of MCMC steps between reruns but a different +##' pattern. ##' ##' @return A `monty_sampler` object, which can be used with ##' [monty_sample] ##' ##' @export -monty_sampler_random_walk <- function(vcv = NULL, boundaries = "reflect") { +monty_sampler_random_walk <- function(vcv = NULL, boundaries = "reflect", + rerun_every = Inf, rerun_random = FALSE) { check_vcv(vcv, allow_3d = TRUE, call = environment()) internal <- new.env() boundaries <- match_value(boundaries, c("reflect", "reject", "ignore")) + + if (!identical(unname(rerun_every), Inf)) { + assert_scalar_positive_integer(rerun_every) + } + assert_scalar_logical(rerun_random) initialise <- function(pars, model, observer, rng) { n_pars <- length(model$parameters) @@ -40,6 +62,10 @@ monty_sampler_random_walk <- function(vcv = NULL, boundaries = "reflect") { ## this is enforced elsewhere stopifnot(model$properties$allow_multiple_parameters) } + + internal$rerun <- + make_rerun(rerun_every, rerun_random, model$properties$is_stochastic) + internal$step <- 0 vcv <- sampler_validate_vcv(vcv, pars) internal$proposal <- @@ -48,6 +74,18 @@ monty_sampler_random_walk <- function(vcv = NULL, boundaries = "reflect") { } step <- function(state, model, observer, rng) { + internal$step <- internal$step + 1 + rerun <- internal$rerun(internal$step, rng) + if (any(rerun)) { + ## This is currently just setup assuming we are not using multiple + ## parameters as currently they cannot be used with stochastic models, + ## while the rerun is only used with stochastic models + state$density <- model$density(state$pars) + if (!is.null(observer)) { + state$observation <- observer$observe(model$model, rng) + } + } + pars_next <- internal$proposal(state$pars, rng) reject_some <- boundaries == "reject" && !all(i <- is_parameters_in_domain(pars_next, model$domain)) @@ -147,3 +185,14 @@ is_parameters_in_domain <- function(x, domain) { if (is.matrix(x)) apply(i, 2, all) else FALSE } } + + +make_rerun <- function(every, random, is_stochastic) { + if (!is_stochastic && !is.finite(every)) { + function(i, rng) rep(FALSE, rng$size()) + } else if (random) { + function(i, rng) rng$random_real(1) < 1 / every + } else { + function(i, rng) rep(i %% every == 0, rng$size()) + } +} diff --git a/man/monty_sampler_random_walk.Rd b/man/monty_sampler_random_walk.Rd index 2b15a835..61f6ad18 100644 --- a/man/monty_sampler_random_walk.Rd +++ b/man/monty_sampler_random_walk.Rd @@ -4,7 +4,12 @@ \alias{monty_sampler_random_walk} \title{Random Walk Sampler} \usage{ -monty_sampler_random_walk(vcv = NULL, boundaries = "reflect") +monty_sampler_random_walk( + vcv = NULL, + boundaries = "reflect", + rerun_every = Inf, + rerun_random = FALSE +) } \arguments{ \item{vcv}{A variance covariance matrix for the proposal.} @@ -23,6 +28,22 @@ the domain. The initial point selected will lie within the domain, as this is enforced by \link{monty_sample}.} + +\item{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 +(\code{Inf}) will never rerun this point, but if you set to 100, then +every 100 steps we run the model on both the proposed \emph{and} previously +accepted parameters before doing the comparison. This may help "unstick" +chains, at the cost of some bias in the results.} + +\item{rerun_random}{Logical, controlling the behaviour of +rerunning (when \code{rerun_every} is finite). The default value of +\code{FALSE} will rerun the model at fixed intervals of iterations +(given by \code{rerun_every}). If \code{TRUE}, then we stochastically +rerun each step with probability of \code{1 / rerun_every}. This gives +the same expected number of MCMC steps between reruns but a different +pattern.} } \value{ A \code{monty_sampler} object, which can be used with From d24b1f45afd2cac966c6f1fcae1e01b6a06ac70a Mon Sep 17 00:00:00 2001 From: edknock Date: Fri, 6 Sep 2024 14:22:29 +0100 Subject: [PATCH 02/14] rerun for nested random walk --- R/sampler-nested-random-walk.R | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/R/sampler-nested-random-walk.R b/R/sampler-nested-random-walk.R index 6f59f907..7506db52 100644 --- a/R/sampler-nested-random-walk.R +++ b/R/sampler-nested-random-walk.R @@ -77,6 +77,11 @@ monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect") { internal <- new.env(parent = emptyenv()) boundaries <- match_value(boundaries, c("reflect", "reject", "ignore")) + + if (!identical(unname(rerun_every), Inf)) { + assert_scalar_positive_integer(rerun_every) + } + assert_scalar_logical(rerun_random) initialise <- function(pars, model, observer, rng) { if (!model$properties$has_parameter_groups) { @@ -88,6 +93,10 @@ monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect") { ## this is enforced elsewhere stopifnot(model$properties$allow_multiple_parameters) } + + internal$rerun <- + make_rerun(rerun_every, rerun_random, model$properties$is_stochastic) + internal$step <- 0 internal$proposal <- nested_proposal(vcv, model$parameter_groups, pars, model$domain, boundaries) @@ -130,6 +139,20 @@ monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect") { ## either changing the behaviour of the step function or swapping in ## a different version. step <- function(state, model, observer, rng) { + internal$step <- internal$step + 1 + rerun <- internal$rerun(internal$step, rng) + if (any(rerun)) { + ## This is currently just setup assuming we are not using multiple + ## parameters as currently they cannot be used with stochastic models, + ## while the rerun is only used with stochastic models + density <- model$density(state$pars, by_group = TRUE) + state$density <- c(density) + internal$density_by_group <- attr(density, "by_group") + if (!is.null(observer)) { + state$observation <- observer$observe(model$model, rng) + } + } + if (!is.null(internal$proposal$base)) { pars_next <- internal$proposal$base(state$pars, rng) From 3ac740d669a77257f272e44dd06317daad604f74 Mon Sep 17 00:00:00 2001 From: edknock Date: Fri, 6 Sep 2024 17:37:45 +0100 Subject: [PATCH 03/14] regenerate utils-assert --- R/import-standalone-utils-assert.R | 50 +++++++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 4 deletions(-) diff --git a/R/import-standalone-utils-assert.R b/R/import-standalone-utils-assert.R index 67adf3bf..22a881e6 100644 --- a/R/import-standalone-utils-assert.R +++ b/R/import-standalone-utils-assert.R @@ -1,6 +1,6 @@ # Standalone file: do not edit by hand -# Source: https://github.com/reside-ic/reside.utils/blob/prototype/R/standalone-utils-assert.R -# Generated by: usethis::use_standalone("reside-ic/reside.utils", "utils-assert", ref = "prototype") +# Source: https://github.com/reside-ic/reside.utils/blob/HEAD/R/standalone-utils-assert.R +# Generated by: usethis::use_standalone("reside-ic/reside.utils", "utils-assert") # ---------------------------------------------------------------------- # # --- @@ -48,7 +48,7 @@ assert_integer <- function(x, name = deparse(substitute(x)), } if (!isTRUE(all.equal(x, rx, tolerance = tolerance))) { cli::cli_abort( - c("Exected '{name}' to be integer", + c("Expected '{name}' to be integer", i = paste("{cli::qty(length(x))}The provided", "{?value was/values were} numeric, but not very close", "to integer values")), @@ -56,7 +56,7 @@ assert_integer <- function(x, name = deparse(substitute(x)), } x <- as.integer(rx) } else { - cli::cli_abort("Exected '{name}' to be integer", call = call, arg = arg) + cli::cli_abort("Expected '{name}' to be integer", call = call, arg = arg) } invisible(x) } @@ -156,6 +156,48 @@ assert_list <- function(x, name = deparse(substitute(x)), arg = name, } +assert_scalar_positive_numeric <- function(x, allow_zero = TRUE, + name = deparse(substitute(x)), + arg = name, call = parent.frame()) { + assert_scalar_numeric(x, name = name, call = call) + if (allow_zero) { + if (x < 0) { + cli::cli_abort("'{name}' must be at least 0", arg = arg, call = call) + } + } else { + if (x <= 0) { + cli::cli_abort("'{name}' must be greater than 0", arg = arg, call = call) + } + } + invisible(x) +} + + +assert_scalar_positive_integer <- function(x, allow_zero = TRUE, + name = deparse(substitute(x)), + tolerance = NULL, arg = name, + call = parent.frame()) { + assert_scalar_integer(x, name, tolerance = tolerance, arg = arg, call = call) + min <- if (allow_zero) 0 else 1 + if (x < min) { + cli::cli_abort("'{name}' must be at least {min}", arg = arg, call = call) + } + invisible(x) +} + + +assert_raw <- function(x, len = NULL, name = deparse(substitute(x)), + arg = name, call = parent.frame()) { + if (!is.raw(x)) { + cli::cli_abort("'{name}' must be a raw vector", arg = arg, call = call) + } + if (!is.null(len)) { + assert_length(x, len, name = name, call = call) + } + invisible(x) +} + + assert_named <- function(x, unique = FALSE, name = deparse(substitute(x)), arg = name, call = parent.frame()) { nms <- names(x) From d90b59f74168cca4acaad2b5c3d3df48245af693 Mon Sep 17 00:00:00 2001 From: edknock Date: Sat, 7 Sep 2024 00:38:47 +0100 Subject: [PATCH 04/14] add rerun tests --- R/sampler-random-walk.R | 2 +- man/monty_sampler_random_walk.Rd | 2 +- tests/testthat/test-sampler-random-walk.R | 16 ++++++++++++++++ 3 files changed, 18 insertions(+), 2 deletions(-) diff --git a/R/sampler-random-walk.R b/R/sampler-random-walk.R index 09310900..8722454d 100644 --- a/R/sampler-random-walk.R +++ b/R/sampler-random-walk.R @@ -26,7 +26,7 @@ ##' @param 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 -##' (`Inf`) will never rerun this point, but if you set to 100, then +##' (`Inf`) will never trigger a rerun, but if you set to 100, then ##' every 100 steps we run the model on both the proposed *and* previously ##' accepted parameters before doing the comparison. This may help "unstick" ##' chains, at the cost of some bias in the results. diff --git a/man/monty_sampler_random_walk.Rd b/man/monty_sampler_random_walk.Rd index 61f6ad18..4d44f68f 100644 --- a/man/monty_sampler_random_walk.Rd +++ b/man/monty_sampler_random_walk.Rd @@ -32,7 +32,7 @@ enforced by \link{monty_sample}.} \item{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 -(\code{Inf}) will never rerun this point, but if you set to 100, then +(\code{Inf}) will never trigger a rerun, but if you set to 100, then every 100 steps we run the model on both the proposed \emph{and} previously accepted parameters before doing the comparison. This may help "unstick" chains, at the cost of some bias in the results.} diff --git a/tests/testthat/test-sampler-random-walk.R b/tests/testthat/test-sampler-random-walk.R index d81c60e1..9fcfde06 100644 --- a/tests/testthat/test-sampler-random-walk.R +++ b/tests/testthat/test-sampler-random-walk.R @@ -342,3 +342,19 @@ test_that("can print a sampler object", { "", fixed = TRUE, all = FALSE) }) + + +test_that("Can rerun a stochastic model", { + set.seed(1) + m <- ex_dust_sir() + 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, + rerun_random = TRUE) + res1 <- monty_sample(m, sampler1, 20) + res2 <- monty_sample(m, sampler2, 20) + + expect_gt(sum(diff(res1$density) != 0), sum(diff(res1$pars[1, , 1]) != 0)) + expect_true(all(diff(res1$density)[seq(1, 20, by = 2)] != 0)) + expect_gt(sum(diff(res2$density) != 0), sum(diff(res2$pars[1, , 1]) != 0)) +}) From 52c279db23dcde7f3a3acc17ba9d5815aeec8072 Mon Sep 17 00:00:00 2001 From: edknock Date: Sat, 7 Sep 2024 00:40:12 +0100 Subject: [PATCH 05/14] bump version number --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index ce7ac01b..a0df548c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.2.3 +Version: 0.2.4 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), From 472561b31f04790c7ebd892c8ef2b2fce7519896 Mon Sep 17 00:00:00 2001 From: edknock Date: Mon, 9 Sep 2024 09:48:02 +0100 Subject: [PATCH 06/14] rerun inputs for nested rw --- R/sampler-nested-random-walk.R | 4 +++- man/monty_sampler_nested_random_walk.Rd | 23 ++++++++++++++++++++++- 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/R/sampler-nested-random-walk.R b/R/sampler-nested-random-walk.R index 7506db52..0e29058e 100644 --- a/R/sampler-nested-random-walk.R +++ b/R/sampler-nested-random-walk.R @@ -49,7 +49,9 @@ ##' [monty_sample] ##' ##' @export -monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect") { +monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect", + rerun_every = Inf, + rerun_random = FALSE) { if (!is.list(vcv)) { cli::cli_abort( "Expected a list for 'vcv'", diff --git a/man/monty_sampler_nested_random_walk.Rd b/man/monty_sampler_nested_random_walk.Rd index e7ceecd7..58642a5f 100644 --- a/man/monty_sampler_nested_random_walk.Rd +++ b/man/monty_sampler_nested_random_walk.Rd @@ -4,7 +4,12 @@ \alias{monty_sampler_nested_random_walk} \title{Nested Random Walk Sampler} \usage{ -monty_sampler_nested_random_walk(vcv, boundaries = "reflect") +monty_sampler_nested_random_walk( + vcv, + boundaries = "reflect", + rerun_every = Inf, + rerun_random = FALSE +) } \arguments{ \item{vcv}{A list of variance covariance matrices. We expect this @@ -25,6 +30,22 @@ the domain. The initial point selected will lie within the domain, as this is enforced by \link{monty_sample}.} + +\item{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 +(\code{Inf}) will never trigger a rerun, but if you set to 100, then +every 100 steps we run the model on both the proposed \emph{and} previously +accepted parameters before doing the comparison. This may help "unstick" +chains, at the cost of some bias in the results.} + +\item{rerun_random}{Logical, controlling the behaviour of +rerunning (when \code{rerun_every} is finite). The default value of +\code{FALSE} will rerun the model at fixed intervals of iterations +(given by \code{rerun_every}). If \code{TRUE}, then we stochastically +rerun each step with probability of \code{1 / rerun_every}. This gives +the same expected number of MCMC steps between reruns but a different +pattern.} } \value{ A \code{monty_sampler} object, which can be used with From 806957c331477e8f5eab41fa50b17f898fad39f9 Mon Sep 17 00:00:00 2001 From: edknock Date: Tue, 10 Sep 2024 11:46:45 +0100 Subject: [PATCH 07/14] begin stochastic nested example --- tests/testthat/helper-monty.R | 42 +++++++++++++++++++ .../test-sampler-nested-random-walk.R | 18 ++++++++ 2 files changed, 60 insertions(+) diff --git a/tests/testthat/helper-monty.R b/tests/testthat/helper-monty.R index 9acd5004..140fa0bc 100644 --- a/tests/testthat/helper-monty.R +++ b/tests/testthat/helper-monty.R @@ -81,6 +81,48 @@ 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") diff --git a/tests/testthat/test-sampler-nested-random-walk.R b/tests/testthat/test-sampler-nested-random-walk.R index 8b22fa9c..0137b828 100644 --- a/tests/testthat/test-sampler-nested-random-walk.R +++ b/tests/testthat/test-sampler-nested-random-walk.R @@ -246,3 +246,21 @@ test_that("can run nested random walk sampler with rejecting boundaries res2 <- monty_sample(m, sampler, 100, n_chains = 3, runner = runner) expect_equal(res1, res2) }) + + +test_that("Can rerun a nested stochastic model", { + set.seed(1) + ng <- 5 + m <- ex_simple_nested_with_base_stochastic(ng) + vcv <- list(base = diag(1), groups = rep(list(diag(1)), ng)) + sampler1 <- monty_sampler_nested_random_walk(vcv, rerun_every = 2) + sampler2 <- monty_sampler_nested_random_walk(vcv, rerun_every = 2, + rerun_random = TRUE) + res1 <- monty_sample(m, sampler1, 20) + res2 <- monty_sample(m, sampler2, 20) + + expect_gt(sum(diff(res1$density) != 0), sum(diff(res1$pars[1, , 1]) != 0)) + expect_true(all(diff(res1$density)[seq(1, 20, by = 2)] != 0)) + expect_gt(sum(diff(res2$density) != 0), sum(diff(res2$pars[1, , 1]) != 0)) +}) + From e3e44ac457954c6b73d1944243663985c6589f1e Mon Sep 17 00:00:00 2001 From: edknock Date: Wed, 2 Oct 2024 13:17:07 +0100 Subject: [PATCH 08/14] bump version number --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9146d764..1a8d44af 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.2.7 +Version: 0.2.8 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), From 15fb107f4b72fe2a54d8ad4b39eaecab42bbd191 Mon Sep 17 00:00:00 2001 From: edknock Date: Wed, 2 Oct 2024 14:21:07 +0100 Subject: [PATCH 09/14] bump version number --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1a8d44af..0cc74a0b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.2.8 +Version: 0.2.9 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), From 67eeacc7328f9fdf2930b7fea646915ea59f0810 Mon Sep 17 00:00:00 2001 From: edknock Date: Wed, 23 Oct 2024 09:14:53 +0100 Subject: [PATCH 10/14] bump version number --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 528079fe..3460b002 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.2.22 +Version: 0.2.23 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), From c78e34c8eae25fbe51f77e32d4fe3873e53c8266 Mon Sep 17 00:00:00 2001 From: edknock Date: Wed, 23 Oct 2024 09:35:23 +0100 Subject: [PATCH 11/14] remove rerun from nested random walk for time being --- R/sampler-nested-random-walk.R | 24 +------------------ .../test-sampler-nested-random-walk.R | 18 -------------- 2 files changed, 1 insertion(+), 41 deletions(-) diff --git a/R/sampler-nested-random-walk.R b/R/sampler-nested-random-walk.R index e8710858..04447f94 100644 --- a/R/sampler-nested-random-walk.R +++ b/R/sampler-nested-random-walk.R @@ -49,9 +49,7 @@ ##' [monty_sample] ##' ##' @export -monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect", - rerun_every = Inf, - rerun_random = FALSE) { +monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect") { if (!is.list(vcv)) { cli::cli_abort( "Expected a list for 'vcv'", @@ -79,11 +77,6 @@ monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect", internal <- new.env(parent = emptyenv()) boundaries <- match_value(boundaries, c("reflect", "reject", "ignore")) - - if (!identical(unname(rerun_every), Inf)) { - assert_scalar_positive_integer(rerun_every) - } - assert_scalar_logical(rerun_random) initialise <- function(pars, model, rng) { if (!model$properties$has_parameter_groups) { @@ -96,10 +89,6 @@ monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect", stopifnot(model$properties$allow_multiple_parameters) } - internal$rerun <- - make_rerun(rerun_every, rerun_random, model$properties$is_stochastic) - internal$step <- 0 - vcv <- sampler_validate_nested_vcv(vcv, model$parameter_groups, pars) internal$proposal <- nested_proposal(vcv, model$parameter_groups, pars, @@ -144,17 +133,6 @@ monty_sampler_nested_random_walk <- function(vcv, boundaries = "reflect", ## either changing the behaviour of the step function or swapping in ## a different version. step <- function(state, model, rng) { - internal$step <- internal$step + 1 - rerun <- internal$rerun(internal$step, rng) - if (any(rerun)) { - ## This is currently just setup assuming we are not using multiple - ## parameters as currently they cannot be used with stochastic models, - ## while the rerun is only used with stochastic models - density <- model$density(state$pars, by_group = TRUE) - state$density <- c(density) - internal$density_by_group <- attr(density, "by_group") - } - if (!is.null(internal$proposal$base)) { pars_next <- internal$proposal$base(state$pars, rng) diff --git a/tests/testthat/test-sampler-nested-random-walk.R b/tests/testthat/test-sampler-nested-random-walk.R index 141e1f90..6ba86111 100644 --- a/tests/testthat/test-sampler-nested-random-walk.R +++ b/tests/testthat/test-sampler-nested-random-walk.R @@ -249,21 +249,3 @@ test_that("can run nested random walk sampler with rejecting boundaries res2 <- monty_sample(m, sampler, 100, n_chains = 3, runner = runner) expect_equal(res1, res2) }) - - -test_that("Can rerun a nested stochastic model", { - set.seed(1) - ng <- 5 - m <- ex_simple_nested_with_base_stochastic(ng) - vcv <- list(base = diag(1), groups = rep(list(diag(1)), ng)) - sampler1 <- monty_sampler_nested_random_walk(vcv, rerun_every = 2) - sampler2 <- monty_sampler_nested_random_walk(vcv, rerun_every = 2, - rerun_random = TRUE) - res1 <- monty_sample(m, sampler1, 20) - res2 <- monty_sample(m, sampler2, 20) - - expect_gt(sum(diff(res1$density) != 0), sum(diff(res1$pars[1, , 1]) != 0)) - expect_true(all(diff(res1$density)[seq(1, 20, by = 2)] != 0)) - expect_gt(sum(diff(res2$density) != 0), sum(diff(res2$pars[1, , 1]) != 0)) -}) - From 5491266517e71ebae3ba82199c379e21a08f0696 Mon Sep 17 00:00:00 2001 From: edknock Date: Wed, 23 Oct 2024 09:37:12 +0100 Subject: [PATCH 12/14] recompile docs --- man/monty_rng.Rd | 53 +++++++++++++++++++++++++ man/monty_sampler_nested_random_walk.Rd | 23 +---------- 2 files changed, 54 insertions(+), 22 deletions(-) diff --git a/man/monty_rng.Rd b/man/monty_rng.Rd index 4642f365..16a16425 100644 --- a/man/monty_rng.Rd +++ b/man/monty_rng.Rd @@ -143,6 +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-beta_binomial_ab}{\code{monty_rng$beta_binomial_ab()}} +\item \href{#method-monty_rng-beta_binomial_prob}{\code{monty_rng$beta_binomial_prob()}} \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()}} @@ -339,6 +341,57 @@ Generate \code{n} numbers from a 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-beta_binomial_ab}{}}} +\subsection{Method \code{beta_binomial_ab()}}{ +Generate \code{n} numbers from a beta-binomial distribution +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{monty_rng$beta_binomial_ab(n, size, a, b, 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 number of trials (zero or more, length 1 or n)} + +\item{\code{a}}{The first shape parameter (zero or more, length 1 or n)} + +\item{\code{b}}{The second shape parameter (zero or more, 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-beta_binomial_prob}{}}} +\subsection{Method \code{beta_binomial_prob()}}{ +Generate \code{n} numbers from a beta-binomial distribution +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{monty_rng$beta_binomial_prob(n, size, prob, rho, 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 number of trials (zero or more, length 1 or n)} + +\item{\code{prob}}{The mean probability of sucess on each trial +(between 0 and 1, length 1 or n)} + +\item{\code{rho}}{The dispersion parameter (between 0 and 1, length 1 or n)} + \item{\code{n_threads}}{Number of threads to use; see Details} } \if{html}{\out{
}} diff --git a/man/monty_sampler_nested_random_walk.Rd b/man/monty_sampler_nested_random_walk.Rd index 58642a5f..e7ceecd7 100644 --- a/man/monty_sampler_nested_random_walk.Rd +++ b/man/monty_sampler_nested_random_walk.Rd @@ -4,12 +4,7 @@ \alias{monty_sampler_nested_random_walk} \title{Nested Random Walk Sampler} \usage{ -monty_sampler_nested_random_walk( - vcv, - boundaries = "reflect", - rerun_every = Inf, - rerun_random = FALSE -) +monty_sampler_nested_random_walk(vcv, boundaries = "reflect") } \arguments{ \item{vcv}{A list of variance covariance matrices. We expect this @@ -30,22 +25,6 @@ the domain. The initial point selected will lie within the domain, as this is enforced by \link{monty_sample}.} - -\item{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 -(\code{Inf}) will never trigger a rerun, but if you set to 100, then -every 100 steps we run the model on both the proposed \emph{and} previously -accepted parameters before doing the comparison. This may help "unstick" -chains, at the cost of some bias in the results.} - -\item{rerun_random}{Logical, controlling the behaviour of -rerunning (when \code{rerun_every} is finite). The default value of -\code{FALSE} will rerun the model at fixed intervals of iterations -(given by \code{rerun_every}). If \code{TRUE}, then we stochastically -rerun each step with probability of \code{1 / rerun_every}. This gives -the same expected number of MCMC steps between reruns but a different -pattern.} } \value{ A \code{monty_sampler} object, which can be used with From bcd1e09528c535dc9b08302e2afaf1c79566ecb2 Mon Sep 17 00:00:00 2001 From: edknock Date: Wed, 23 Oct 2024 10:18:13 +0100 Subject: [PATCH 13/14] make rerun_random = TRUE default --- R/sampler-random-walk.R | 12 ++++++------ man/monty_sampler_random_walk.Rd | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/R/sampler-random-walk.R b/R/sampler-random-walk.R index d338e440..e66605f3 100644 --- a/R/sampler-random-walk.R +++ b/R/sampler-random-walk.R @@ -32,11 +32,11 @@ ##' chains, at the cost of some bias in the results. ##' ##' @param rerun_random Logical, controlling the behaviour of -##' rerunning (when `rerun_every` is finite). The default value of -##' `FALSE` will rerun the model at fixed intervals of iterations -##' (given by `rerun_every`). If `TRUE`, then we stochastically -##' rerun each step with probability of `1 / rerun_every`. This gives -##' the same expected number of MCMC steps between reruns but a different +##' rerunning (when `rerun_every` is finite). With the default value of +##' `TRUE`, we stochastically rerun at each step with probability of +##' `1 / rerun_every`. If `FALSE` we rerun the model at fixed intervals of +##' iterations (given by `rerun_every`). The two methods give the same +##' expected number of MCMC steps between reruns but a different ##' pattern. ##' ##' @return A `monty_sampler` object, which can be used with @@ -44,7 +44,7 @@ ##' ##' @export monty_sampler_random_walk <- function(vcv = NULL, boundaries = "reflect", - rerun_every = Inf, rerun_random = FALSE) { + rerun_every = Inf, rerun_random = TRUE) { check_vcv(vcv, allow_3d = TRUE, call = environment()) internal <- new.env() diff --git a/man/monty_sampler_random_walk.Rd b/man/monty_sampler_random_walk.Rd index 4d44f68f..d592f9ea 100644 --- a/man/monty_sampler_random_walk.Rd +++ b/man/monty_sampler_random_walk.Rd @@ -8,7 +8,7 @@ monty_sampler_random_walk( vcv = NULL, boundaries = "reflect", rerun_every = Inf, - rerun_random = FALSE + rerun_random = TRUE ) } \arguments{ @@ -38,11 +38,11 @@ accepted parameters before doing the comparison. This may help "unstick" chains, at the cost of some bias in the results.} \item{rerun_random}{Logical, controlling the behaviour of -rerunning (when \code{rerun_every} is finite). The default value of -\code{FALSE} will rerun the model at fixed intervals of iterations -(given by \code{rerun_every}). If \code{TRUE}, then we stochastically -rerun each step with probability of \code{1 / rerun_every}. This gives -the same expected number of MCMC steps between reruns but a different +rerunning (when \code{rerun_every} is finite). With the default value of +\code{TRUE}, we stochastically rerun at each step with probability of +\code{1 / rerun_every}. If \code{FALSE} we rerun the model at fixed intervals of +iterations (given by \code{rerun_every}). The two methods give the same +expected number of MCMC steps between reruns but a different pattern.} } \value{ From a0cd09fe66731ecb000030c1999b8a20004b959e Mon Sep 17 00:00:00 2001 From: edknock Date: Wed, 23 Oct 2024 10:47:36 +0100 Subject: [PATCH 14/14] fix test --- tests/testthat/test-sampler-random-walk.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-sampler-random-walk.R b/tests/testthat/test-sampler-random-walk.R index 8229b709..d789d7b9 100644 --- a/tests/testthat/test-sampler-random-walk.R +++ b/tests/testthat/test-sampler-random-walk.R @@ -330,9 +330,9 @@ test_that("Can rerun a stochastic model", { set.seed(1) 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, - rerun_random = TRUE) + sampler1 <- monty_sampler_random_walk(vcv = vcv, rerun_every = 2, + rerun_random = FALSE) + sampler2 <- monty_sampler_random_walk(vcv = vcv, rerun_every = 2) res1 <- monty_sample(m, sampler1, 20) res2 <- monty_sample(m, sampler2, 20)