version 2.0

DESCRIPTION

@@ -0,0 +1,38 @@
+Package: mcmcensemble
+Title: Ensemble Sampler for Affine-Invariant MCMC
+Version: 2.0
+Authors@R: 
+    c(person(given = "Hugo",
+             family = "Gruson",
+             role = c("cre", "aut", "cph"),
+             email = "hugo.gruson+R@normalesup.org",
+             comment = c(ORCID = "0000-0002-4094-1476")),
+      person(given = "Sanda",
+             family = "Dejanic",
+             role = c("aut", "cph")),
+      person(given = "Andreas",
+             family = "Scheidegger",
+             role = c("aut", "cph"),
+             comment = c(ORCID = "0000-0003-2575-2172")))
+Description: Provides ensemble samplers for
+    affine-invariant Monte Carlo Markov Chain, which allow a faster
+    convergence for badly scaled estimation problems. Two samplers are
+    proposed: the 'differential.evolution' sampler from ter Braak and
+    Vrugt (2008) <doi:10.1007/s11222-008-9104-9> and the 'stretch' sampler
+    from Goodman and Weare (2010) <doi:10.2140/camcos.2010.5.65>.
+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
+Encoding: UTF-8
+RoxygenNote: 7.1.1
+NeedsCompilation: no
+Packaged: 2020-10-09 14:08:44 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

MD5

@@ -0,0 +1,17 @@
+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
+7d50348b7238627891f0c28249748996 *man/s.m.mcmc.Rd
+1495e1617d4670efb997b4dc2f4717f2 *tests/testthat.R
+69eea9edded2570e2630d38a80c2d993 *tests/testthat/test-MCMCEnsemble.R

NAMESPACE

@@ -0,0 +1,6 @@
+# Generated by roxygen2: do not edit by hand
+
+export(MCMCEnsemble)
+export(d.e.mcmc)
+export(s.m.mcmc)
+importFrom(stats,runif)

NEWS.md

@@ -0,0 +1,17 @@
+# mcmcensemble 2.0
+
+## Breaking changes
+
+* The argument names and order in `d.e.mcmc()` and `s.m.mcmc()` now match those
+of `MCMCEnsemble()`
+
+## Other user-facing changes
+
+* coda package is now only in `Suggests`, instead of being a hard dependency
+
+## Dev changes
+
+* this package is now named mcmcensemble
+* roxygen2 documentation now uses markdown syntax
+* this package now has unit and regression tests
+* various parts of the code have been optimized for speed

R/MCMCEnsemble.R

@@ -0,0 +1,127 @@
+#' MCMC ensemble sampler
+#'
+#' Ensemble Markov Chain Monte Carlo sampler with different strategies to
+#' generate proposals. Either the *stretch move* as proposed by Goodman and
+#' Weare (2010), or a *differential evolution jump move* similar to Braak and
+#' Vrugt (2008).
+#'
+#' @param f function that returns a single scalar value proportional to the log
+#'   probability density to sample from.
+#' @param lower.inits vector specifying for each parameter the lower value the
+#'   initial distribution.
+#' @param upper.inits vector specifying for each parameter the upper value the
+#'   initial distribution.
+#' @param max.iter maximum number of function evaluations
+#' @param n.walkers number of walkers (ensemble size)
+#' @param method method for proposal generation, either `"stretch"`, or
+#'   `"differential.evolution"`.
+#' @param coda logical. Should the samples be returned as [coda::mcmc.list]
+#'   object? (defaults to `FALSE`)
+#' @param ... further arguments passed to `f`
+#'
+#' @return
+#' * if `coda = FALSE` a list with:
+#' - *samples*: A three dimensional array of samples with dimensions `walker` x
+#' `generation` x `parameter`
+#' - *log.p*: A matrix with the log density evaluate for each walker at each
+#' generation.
+#' * if `coda = TRUE` a list with:
+#' - *samples*: A object of class [coda::mcmc.list] containing all samples.
+#' - *log.p*: A matrix with the log density evaluate for each walker at each
+#' generation.
+#'
+#' @examples
+#' ## a log-pdf to sample from
+#' p.log <- function(x) {
+#'     B <- 0.03                              # controls 'bananacity'
+#'     -x[1]^2/200 - 1/2*(x[2]+B*x[1]^2-100*B)^2
+#' }
+#'
+#' ## 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")
+#' 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",
+#'                      coda=TRUE)
+#'
+#' summary(res2$samples)
+#' plot(res2$samples)
+#'
+#'
+#' ## 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,
+#'                      method="differential.evolution", coda=TRUE)
+#'
+#' summary(res3$samples)
+#' plot(res3$samples)
+#'
+#' @export
+#'
+#' @references
+#' - 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.
+#' Communications in Applied Mathematics and Computational Science, 5(1), 65–80,
+#' \doi{10.2140/camcos.2010.5.65}
+#'
+MCMCEnsemble <- function(f, lower.inits, upper.inits,
+                         max.iter, n.walkers = 10 * length(lower.inits),
+                         method = c("stretch", "differential.evolution"),
+                         coda = FALSE, ...) {
+  if (length(lower.inits) != length(upper.inits)) {
+    stop("The length of 'lower.inits' and 'lower.inits' is must be identical!")
+  }
+
+  n.dim <- length(lower.inits)
+  init.range <- cbind(lower.inits, upper.inits)
+
+  ## run mcmc
+  method <- match.arg(method)
+  message("Using ", method, " move with ", n.walkers, " walkers.")
+
+  if (method == "differential.evolution") {
+    res <- d.e.mcmc(f, lower.inits, upper.inits, max.iter, n.walkers, ...)
+  }
+  if (method == "stretch") {
+    res <- s.m.mcmc(f, lower.inits, upper.inits, max.iter, n.walkers, ...)
+  }
+
+  ## add names
+  if (is.null(names(lower.inits))) {
+    pnames <- paste0("para_", 1:n.dim)
+  } else {
+    pnames <- names(lower.inits)
+  }
+  dimnames(res$samples) <- list(
+    paste0("walker_", 1:n.walkers),
+    paste0("generation_", 1:dim(res$samples)[2]),
+    pnames
+  )
+
+  dimnames(res$log.p) <- list(
+    paste0("walker_", 1:n.walkers),
+    paste0("generation_", 1:dim(res$samples)[2])
+  )
+
+  ## convert to coda object
+  if (coda) {
+
+    if (!requireNamespace("coda", quietly = TRUE)) {
+      stop("Package \"coda\" needed for projection plots. Please install it.",
+           call. = FALSE
+      )
+    }
+
+    ll <- lapply(seq_len(n.walkers),
+                 function(w) coda::as.mcmc(res$samples[w, , ]))
+    res <- list(samples = coda::as.mcmc.list(ll), log.p = res$log.p)
+  }
+
+  res
+}

R/MCMCEnsembleSampler-package.R

@@ -0,0 +1,8 @@
+#' MCMCEnsembleSampler
+#'
+#' This package implements a particle Monte Carlo Markov Chain sampler
+#' with two different ways of creating proposals.
+#'
+#' @keywords internal
+"_PACKAGE"
+

R/d.e.mcmc.R

@@ -0,0 +1,88 @@
+#' 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)
+#'
+#' @inheritParams MCMCEnsemble
+#'
+#' @author Sanda Dejanic
+#'
+#' @return List containing:
+#' - `samples[n.walkers,chain.length,n.dim]`
+#' - `log.p[n.walkers,chain.length]`
+#'
+#' @importFrom stats runif
+#'
+#' @export
+#'
+#' @references
+#' 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}
+#' .
+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))
+
+  ensemble.old[1, ] <- runif(n.dim, lower.inits, upper.inits)
+  logres <- f(ensemble.old[1, ], ...)
+  if (length(logres) != 1 || !is.numeric(logres)) {
+    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[, 1] <- log.p.old
+  samples[, 1, ] <- ensemble.old
+
+  ## 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]
+    }
+  }
+
+  mcmc.list <- list(samples = samples, log.p = log.p)
+
+  return(mcmc.list)
+}

R/s.m.mcmc.R

@@ -0,0 +1,79 @@
+#' MCMC Ensemble sampler with the stretch move (emcee)
+#'
+#' Markov Chain Monte Carlo sampler: using the stretch move (implementation of
+#' the Goodman and Ware emcee)
+#'
+#' @inheritParams MCMCEnsemble
+#'
+#' @author Sanda Dejanic
+#'
+#' @inherit d.e.mcmc return
+#'
+#' @importFrom stats runif
+#'
+#' @export
+#'
+#' @references
+#' 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}
+#'
+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))
+
+  ensemble.old[1, ] <- runif(n.dim, lower.inits, upper.inits)
+  logres <- f(ensemble.old[1, ], ...)
+  if (length(logres) != 1 || !is.numeric(logres)) {
+    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[, 1] <- log.p.old
+  samples[, 1, ] <- ensemble.old
+
+  ## 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]
+    }
+  }
+
+  mcmc.list <- list(samples = samples, log.p = log.p)
+
+  return(mcmc.list)
+}