mcmcensemble
This R package 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) and the ‘stretch’ sampler from Goodman and Weare (2010).
Installation
# install.packages("remotes")
remotes::install_github("Bisaloo/mcmcensemble")
Usage
library(mcmcensemble)
## 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")
#> Using stretch move with 10 walkers.
str(res1)
#> List of 2
#> $ samples: num [1:10, 1:300, 1:2] 0.14 0.665 0.995 0.653 0.476 ...
#> ..- attr(*, "dimnames")=List of 3
#> .. ..$ : chr [1:10] "walker_1" "walker_2" "walker_3" "walker_4" ...
#> .. ..$ : chr [1:300] "generation_1" "generation_2" "generation_3" "generation_4" ...
#> .. ..$ : chr [1:2] "a" "b"
#> $ log.p : num [1:10, 1:300] -3.59 -3.48 -2.41 -3.44 -2.44 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:10] "walker_1" "walker_2" "walker_3" "walker_4" ...
#> .. ..$ : chr [1:300] "generation_1" "generation_2" "generation_3" "generation_4" ...
If the coda package is
installed, you can then use the coda = TRUE
argument to get objects of
class mcmc.list
. The coda package then allows you to call summary()
and plot()
to get informative and nicely formatted results and plots:
## 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)
#> Using stretch move with 10 walkers.
summary(res2$samples)
#>
#> Iterations = 1:300
#> Thinning interval = 1
#> Number of chains = 10
#> Sample size per chain = 300
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> a -0.99617 9.404 0.17169 1.1082
#> b 0.08097 3.802 0.06941 0.5324
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> a -20.18 -6.6065 0.01745 5.148 17.65
#> b -10.22 -0.9395 1.25350 2.573 4.29
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)
#> Using differential.evolution move with 10 walkers.
summary(res3$samples)
#>
#> Iterations = 1:300
#> Thinning interval = 1
#> Number of chains = 10
#> Sample size per chain = 300
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> a -0.6493 8.119 0.14824 0.631
#> b 1.0218 2.408 0.04397 0.218
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> a -14.947 -6.59386 -1.001 5.112 16.622
#> b -5.064 -0.08247 1.645 2.747 4.148
plot(res3$samples)
Similar projects
The methods used in this package also have (independent) implementations in other languages:
- emcee in Python (https://doi.org/10.21105/joss.01864)
- gwmcmc in Matlab (https://github.com/grinsted/gwmcmc)