version 2.1

cran · Jan 7, 2021 · 9a52c22 · 9a52c22
1 parent 6e4ec6a
commit 9a52c22
Show file tree

Hide file tree

Showing 24 changed files with 2,402 additions and 1,852 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: mcmcensemble
 Title: Ensemble Sampler for Affine-Invariant MCMC
-Version: 2.0
+Version: 2.1
 Authors@R: 
     c(person(given = "Hugo",
              family = "Gruson",
@@ -24,15 +24,18 @@ License: GPL-2
 URL: https://github.com/Bisaloo/mcmcensemble,
         https://bisaloo.github.io/mcmcensemble/
 BugReports: https://github.com/Bisaloo/mcmcensemble/issues
-Suggests: coda, mockery, testthat
+Imports: future.apply, progressr
+Suggests: coda, mockery, testthat (>= 3.0.0), knitr, rmarkdown
 Encoding: UTF-8
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.1.9000
+Config/testthat/edition: 3
+VignetteBuilder: knitr
 NeedsCompilation: no
-Packaged: 2020-10-09 14:08:44 UTC; hugo
+Packaged: 2021-01-07 10:28:21 UTC; hugo
 Author: Hugo Gruson [cre, aut, cph] (<https://orcid.org/0000-0002-4094-1476>),
   Sanda Dejanic [aut, cph],
   Andreas Scheidegger [aut, cph]
     (<https://orcid.org/0000-0003-2575-2172>)
 Maintainer: Hugo Gruson <hugo.gruson+R@normalesup.org>
 Repository: CRAN
-Date/Publication: 2020-10-17 12:20:02 UTC
+Date/Publication: 2021-01-07 16:10:13 UTC
diff --git a/MD5 b/MD5
@@ -1,17 +1,23 @@
-341c7b0b77949d6b331971a592eaec69 *DESCRIPTION
-23ccf03953716b727e4eec95b4287cca *NAMESPACE
-6dffced3e7e90013ef7e0fa117714635 *NEWS.md
-28de20233efa24c56f88d56b7fcf8516 *R/MCMCEnsemble.R
-5a2436fe6a4782d4b94f91bf9ee820ce *R/MCMCEnsembleSampler-package.R
-83b75a5b140fbb9ae445c3d9c7d92dac *R/d.e.mcmc.R
-545542fbe350f0807f1569f70c4bf2f9 *R/s.m.mcmc.R
-33bdc1720795eee7735a29e32346e074 *README.md
-d5a1a0aa189c8f0e3f033c3c8c7f8a3e *build/partial.rdb
-f9793199ca888f22c690a1ea133ea095 *man/MCMCEnsemble.Rd
-0724f101c1699f085ae6d911db29bb41 *man/d.e.mcmc.Rd
-4da670a6b2d943e54ed52150dc06790a *man/figures/README-unnamed-chunk-4-1.svg
-e85d4784dcd3b3b339cc30ece62bee94 *man/figures/README-unnamed-chunk-4-2.svg
-b24875fcd70e92ecb34171d95e541771 *man/mcmcensemble-package.Rd
+edf684cb7a900d99f175181bbc87abf9 *DESCRIPTION
+f3eaa8908ea4c90dfdeeaaa1fe8b7ada *NAMESPACE
+d33c20848e072bec4a620ae8b224736c *NEWS.md
+6238383143ab613ece4f9d211d59b2f1 *R/MCMCEnsemble.R
+f6d18adfdbf92714f2e96dab29967ff5 *R/d.e.mcmc.R
+5a2436fe6a4782d4b94f91bf9ee820ce *R/mcmcensemble-package.R
+d433b871b1d7812d58f177eefda90441 *R/s.m.mcmc.R
+2f9260bbe6f9092c2e814d66d212d934 *README.md
+33b41c2da5f1cbced0be5e34e6cc19c5 *build/partial.rdb
+14406e0e924f7e5d059bd46583fb664d *build/vignette.rds
+cc80efe04615c251caae9c85ce89cfd1 *inst/doc/faq.R
+0d4312a5f45639b1d3d7f472608d2cd8 *inst/doc/faq.Rmd
+87cb8103e9d3d5d572a24589549a8bf6 *inst/doc/faq.html
+18b8b8e1a04fd0cb0dcdd62586e4bdd8 *man/MCMCEnsemble.Rd
+38c970375477a873d8ac745129091eaa *man/d.e.mcmc.Rd
+61ad5e732f3463fbd9f4f7ffd139d222 *man/figures/README-example-de-1.svg
+79a1d4e6e0cb281a69cd43156edfbac8 *man/figures/README-example-stretch-1.svg
+23c04d0dcaf80725c567cfb046b801ab *man/mcmcensemble-package.Rd
 7d50348b7238627891f0c28249748996 *man/s.m.mcmc.Rd
 1495e1617d4670efb997b4dc2f4717f2 *tests/testthat.R
-69eea9edded2570e2630d38a80c2d993 *tests/testthat/test-MCMCEnsemble.R
+238c0972382e28e608684e8513245b81 *tests/testthat/test-MCMCEnsemble.R
+c8933ae3a93c2fa4e1281dc42d1bd552 *tests/testthat/test-convergence.R
+0d4312a5f45639b1d3d7f472608d2cd8 *vignettes/faq.Rmd
diff --git a/NAMESPACE b/NAMESPACE
@@ -3,4 +3,6 @@
 export(MCMCEnsemble)
 export(d.e.mcmc)
 export(s.m.mcmc)
+importFrom(future.apply,future_apply)
+importFrom(progressr,progressor)
 importFrom(stats,runif)
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,22 @@
+# mcmcensemble 2.1
+
+## Major changes
+
+* the ensemble sampling can now be parallelized with the future framework. Check
+the [README](https://bisaloo.github.io/mcmcensemble/) for more information
+
+## Other user-facing changes
+
+* very large log.p differences between chains do not cause them to be
+stuck any more
+* addition of a new vignette listing frequently asked questions (with their
+answer)
+
+## Dev changes
+
+* new test to make sure the chains converge as expected
+* performance improvements
+
 # mcmcensemble 2.0
 
 ## Breaking changes

diff --git a/R/MCMCEnsemble.R b/R/MCMCEnsemble.R
@@ -39,13 +39,13 @@
 #'
 #' ## use stretch move
 #' res1 <- MCMCEnsemble(p.log, lower.inits=c(a=0, b=0), upper.inits=c(a=1, b=1),
-#'                      max.iter=3000, n.walkers=10, method="stretch")
+#'                      max.iter=300, n.walkers=10, method="stretch")
 #' str(res1)
 #'
 #'
 #' ## use stretch move, return samples as 'coda' object
 #' res2 <- MCMCEnsemble(p.log, lower.inits=c(a=0, b=0), upper.inits=c(a=1, b=1),
-#'                      max.iter=3000, n.walkers=10, method="stretch",
+#'                      max.iter=300, n.walkers=10, method="stretch",
 #'                      coda=TRUE)
 #'
 #' summary(res2$samples)
@@ -54,7 +54,7 @@
 #'
 #' ## use different evolution move, return samples as 'coda' object
 #' res3 <- MCMCEnsemble(p.log, lower.inits=c(a=0, b=0), upper.inits=c(a=1, b=1),
-#'                      max.iter=3000, n.walkers=10,
+#'                      max.iter=300, n.walkers=10,
 #'                      method="differential.evolution", coda=TRUE)
 #'
 #' summary(res3$samples)
@@ -66,7 +66,7 @@
 #' - ter Braak, C. J. F. and Vrugt, J. A. (2008) Differential Evolution Markov
 #' Chain with snooker updater and fewer chains. Statistics and Computing, 18(4),
 #' 435–446, \doi{10.1007/s11222-008-9104-9}
-#' -  Goodman, J. and Weare, J. (2010) Ensemble samplers with affine invariance.
+#' - Goodman, J. and Weare, J. (2010) Ensemble samplers with affine invariance.
 #' Communications in Applied Mathematics and Computational Science, 5(1), 65–80,
 #' \doi{10.2140/camcos.2010.5.65}
 #'
@@ -113,8 +113,10 @@ MCMCEnsemble <- function(f, lower.inits, upper.inits,
   if (coda) {
 
     if (!requireNamespace("coda", quietly = TRUE)) {
-      stop("Package \"coda\" needed for projection plots. Please install it.",
-           call. = FALSE
+      stop(
+        "Package 'coda' needed for to create coda objects. ",
+        "Please install it or use `coda = TRUE`.",
+        call. = FALSE
       )
     }
 

diff --git a/R/d.e.mcmc.R b/R/d.e.mcmc.R
@@ -1,7 +1,7 @@
 #' MCMC Ensemble sampler with the differential evolution jump move
 #'
 #' Markov Chain Monte Carlo sampler: using the differential evolution jump move
-#' (implementation of the Ter Braak differential evolution)
+#' (implementation of the ter Braak differential evolution)
 #'
 #' @inheritParams MCMCEnsemble
 #'
@@ -19,70 +19,78 @@
 #' ter Braak, C. J. F. and Vrugt, J. A. (2008) Differential Evolution Markov
 #' Chain with snooker updater and fewer chains. Statistics and Computing,
 #' 18(4), 435–446, \doi{10.1007/s11222-008-9104-9}
-#' .
+#'
+#' @importFrom future.apply future_apply
+#'
 d.e.mcmc <- function(f, lower.inits, upper.inits, max.iter, n.walkers, ...) {
 
   n.dim <- length(lower.inits)
   ## initial values
 
   chain.length <- max.iter %/% n.walkers
 
-  log.p <- matrix(NA_real_, nrow = n.walkers, ncol = chain.length)
-  log.p.old <- rep(NA_real_, n.walkers)
-  ensemble.old <- matrix(NA_real_, nrow = n.walkers, ncol = n.dim)
-  ensemble.new <- matrix(NA_real_, nrow = n.walkers, ncol = n.dim)
-  samples <- array(NA_real_, dim = c(n.walkers, chain.length, n.dim))
+  p <- progressor(chain.length)
+
+  ensemble.old <- matrix(
+    runif(n.dim*n.walkers, lower.inits, upper.inits),
+    nrow = n.walkers,
+    ncol = n.dim,
+    byrow = TRUE
+  )
+
+  log.p.old <- future_apply(ensemble.old, 1, f, ..., future.seed = TRUE)
 
-  ensemble.old[1, ] <- runif(n.dim, lower.inits, upper.inits)
-  logres <- f(ensemble.old[1, ], ...)
-  if (length(logres) != 1 || !is.numeric(logres)) {
+  if (!is.vector(log.p.old, mode = "numeric")) {
     stop("Function 'f' should return a numeric of length 1", call. = FALSE)
   }
-  log.p.old[1] <- logres
 
-  for (k in 2:n.walkers) {
-    ensemble.old[k, ] <- runif(n.dim, lower.inits, upper.inits)
-    log.p.old[k] <- f(ensemble.old[k, ], ...)
-  }
+  log.p <- matrix(NA_real_, nrow = n.walkers, ncol = chain.length)
+  samples <- array(NA_real_, dim = c(n.walkers, chain.length, n.dim))
 
   log.p[, 1] <- log.p.old
   samples[, 1, ] <- ensemble.old
 
+  p()
+
   ## the loop
 
   for (l in 2:chain.length) {
-    for (n in 1:n.walkers) {
-      z <- 2.38 / sqrt(2 * n.dim)
-      if (l %% 10 == 0) {
-        z <- 1
-      }
-
-      a <- sample((1:n.walkers)[-n], 1)
-      b <- sample((1:n.walkers)[-c(n, a)], 1)
-
-      par.active.1 <- ensemble.old[a, ]
-      par.active.2 <- ensemble.old[b, ]
-
-      ensemble.new[n, ] <- ensemble.old[n, ] + z * (par.active.1 - par.active.2)
-
-      log.p.new <- f(ensemble.new[n, ], ...)
-      if (!is.finite(log.p.new)) {
-        acc <- 0
-      }
-      else {
-        acc <- exp(log.p.new - log.p.old[n])
-      }
-
-      if (acc > runif(1)) {
-        ensemble.old[n, ] <- ensemble.new[n, ]
-        log.p.old[n] <- log.p.new
-      }
-      samples[n, l, ] <- ensemble.old[n, ]
-      log.p[n, l] <- log.p.old[n]
+
+    z <- 2.38 / sqrt(2 * n.dim)
+    if (l %% 10 == 0) {
+      z <- 1
     }
+
+    s <- vapply(
+      seq_len(n.walkers),
+      function(n) sample((1:n.walkers)[-n], 2, replace = FALSE),
+      integer(2)
+    )
+
+    par.active.1 <- ensemble.old[s[1, ], ]
+    par.active.2 <- ensemble.old[s[2, ], ]
+
+    ensemble.new <- ensemble.old + z * (par.active.1 - par.active.2)
+
+    log.p.new <- future_apply(ensemble.new, 1, f, ..., future.seed = TRUE)
+
+    val <- exp(log.p.new - log.p.old)
+
+    # We don't want to get rid of Inf values since +Inf is a valid value to
+    # accept a change. If we forbid, we are actually forbidding large log.p
+    # changes.
+    acc <- !is.na(val) & val > runif(n.walkers)
+
+    ensemble.old[acc, ] <- ensemble.new[acc, ]
+    log.p.old[acc] <- log.p.new[acc]
+
+    samples[, l, ] <- ensemble.old
+    log.p[, l] <- log.p.old
+
+    p()
+
   }
 
-  mcmc.list <- list(samples = samples, log.p = log.p)
+  return(list(samples = samples, log.p = log.p))
 
-  return(mcmc.list)
 }
diff --git a/R/MCMCEnsembleSampler-package.R → R/mcmcensemble-package.R b/R/MCMCEnsembleSampler-package.R → R/mcmcensemble-package.R
diff --git a/R/s.m.mcmc.R b/R/s.m.mcmc.R
@@ -18,62 +18,73 @@
 #' Communications in Applied Mathematics and Computational Science, 5(1), 65–80,
 #' \doi{10.2140/camcos.2010.5.65}
 #'
+#' @importFrom future.apply future_apply
+#' @importFrom progressr progressor
+#'
 s.m.mcmc <- function(f, lower.inits, upper.inits, max.iter, n.walkers, ...) {
 
   n.dim <- length(lower.inits)
   ## initial values
 
   chain.length <- max.iter %/% n.walkers
 
-  log.p <- matrix(NA_real_, nrow = n.walkers, ncol = chain.length)
-  log.p.old <- rep(NA_real_, n.walkers)
-  ensemble.old <- matrix(NA_real_, nrow = n.walkers, ncol = n.dim)
-  ensemble.new <- matrix(NA_real_, nrow = n.walkers, ncol = n.dim)
-  samples <- array(NA_real_, dim = c(n.walkers, chain.length, n.dim))
+  p <- progressor(chain.length)
+
+  ensemble.old <- matrix(
+    runif(n.dim*n.walkers, lower.inits, upper.inits),
+    nrow = n.walkers,
+    ncol = n.dim,
+    byrow = TRUE
+  )
+
+  log.p.old <- future_apply(ensemble.old, 1, f, ..., future.seed = TRUE)
 
-  ensemble.old[1, ] <- runif(n.dim, lower.inits, upper.inits)
-  logres <- f(ensemble.old[1, ], ...)
-  if (length(logres) != 1 || !is.numeric(logres)) {
+  if (!is.vector(log.p.old, mode = "numeric")) {
     stop("Function 'f' should return a numeric of length 1", call. = FALSE)
   }
-  log.p.old[1] <- logres
 
-  for (k in 2:n.walkers) {
-    ensemble.old[k, ] <- runif(n.dim, lower.inits, upper.inits)
-    log.p.old[k] <- f(ensemble.old[k, ], ...)
-  }
+  log.p <- matrix(NA_real_, nrow = n.walkers, ncol = chain.length)
+  samples <- array(NA_real_, dim = c(n.walkers, chain.length, n.dim))
 
   log.p[, 1] <- log.p.old
   samples[, 1, ] <- ensemble.old
 
+  p()
+
   ## the loop
 
   for (l in 2:chain.length) {
-    for (n in 1:n.walkers) {
-      z <- ((runif(1) + 1)^2) / 2
-      a <- sample((1:n.walkers)[-n], 1)
-      par.active <- ensemble.old[a, ]
-
-      ensemble.new[n, ] <- par.active + z * (ensemble.old[n, ] - par.active)
-
-      log.p.new <- f(ensemble.new[n, ], ...)
-      if (!is.finite(log.p.new)) {
-        acc <- 0
-      }
-      else {
-        acc <- z^(n.dim - 1) * exp(log.p.new - log.p.old[n])
-      }
-
-      if (acc > runif(1)) {
-        ensemble.old[n, ] <- ensemble.new[n, ]
-        log.p.old[n] <- log.p.new
-      }
-      samples[n, l, ] <- ensemble.old[n, ]
-      log.p[n, l] <- log.p.old[n]
-    }
+
+    z <- ((runif(n.walkers) + 1)^2) / 2
+
+    a <- vapply(
+      seq_len(n.walkers),
+      function(n) sample((1:n.walkers)[-n], 1),
+      integer(1)
+    )
+
+    par.active <- ensemble.old[a, ]
+
+    ensemble.new <- par.active + z * (ensemble.old - par.active)
+
+    log.p.new <- future_apply(ensemble.new, 1, f, ..., future.seed = TRUE)
+
+    val <- z^(n.dim - 1) * exp(log.p.new - log.p.old)
+
+    # We don't want to get rid of Inf values since +Inf is a valid value to
+    # accept a change. If we forbid, we are actually forbidding large log.p
+    # changes.
+    acc <- !is.na(val) & val > runif(n.walkers)
+
+    ensemble.old[acc, ] <- ensemble.new[acc, ]
+    log.p.old[acc] <- log.p.new[acc]
+
+    samples[, l, ] <- ensemble.old
+    log.p[, l] <- log.p.old
+
+    p()
   }
 
-  mcmc.list <- list(samples = samples, log.p = log.p)
+  return(list(samples = samples, log.p = log.p))
 
-  return(mcmc.list)
 }