Merge pull request #68 from mrc-ide/implement-rerun

Implement rerun for random walk samplers

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 = "[email protected]"),
              person("Wes", "Hinsley", role = "aut"),

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 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.
+##'
+##' @param rerun_random Logical, controlling the behaviour of
+##'   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
 ##'   [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 = TRUE) {
   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, 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,15 @@ monty_sampler_random_walk <- function(vcv = NULL, boundaries = "reflect") {
   }
 
   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
+      state$density <- model$density(state$pars)
+    }
+
     pars_next <- internal$proposal(state$pars, rng)
     reject_some <- boundaries == "reject" &&
       !all(i <- is_parameters_in_domain(pars_next, model$domain))
@@ -147,3 +182,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())
+  }
+}

tests/testthat/test-sampler-random-walk.R

@@ -324,3 +324,19 @@ test_that("can print a sampler object", {
     "<monty_sampler: Random walk (monty_random_walk)>",
     fixed = TRUE, all = FALSE)
 })
+
+
+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,
+                                        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)
+
+  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))
+})