Frequently Asked Questions

library(mcmcensemble)

Can estimations go beyond lower.inits and upper.inits?

YES

## 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
}

set.seed(20201209)

res1 <- MCMCEnsemble(
  p.log,
  lower.inits = c(a = -10, b = -10), upper.inits = c(a = -5, b = -5),
  max.iter = 3000, n.walkers = 10,
  method = "stretch",
  coda = TRUE
)
#> Using stretch move with 10 walkers.

summary(res1$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 -3.4972 9.084  0.16584         1.4141
#> b -0.1255 4.116  0.07515         0.5463
#> 
#> 2. Quantiles for each variable:
#> 
#>      2.5%     25%    50%   75%  97.5%
#> a -20.642 -10.134 -2.887 3.098 12.403
#> b  -9.974  -1.766  1.148 2.675  4.401
plot(res1$samples)
res2 <- MCMCEnsemble(
  p.log,
  lower.inits = c(a = -10, b = -10), upper.inits = c(a = -5, b = -5),
  max.iter = 3000, n.walkers = 10,
  method = "differential.evolution",
  coda = TRUE
)
#> Using differential.evolution 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.14194 9.874  0.18028         0.9203
#> b -0.01001 3.697  0.06751         0.4142
#> 
#> 2. Quantiles for each variable:
#> 
#>      2.5%    25%     50%   75% 97.5%
#> a -19.189 -6.265 -0.1711 7.735 18.07
#> b  -8.556 -1.963  1.2788 2.687  4.34
plot(res2$samples)

How to restrict the possible parameter range?

There is no built-in way to define hard limits for the parameter and make sure they never go outside of this range.

The recommended way to address this issue is to handle these cases in the function f you provide.

For example, to keep parameters in the 0-1 range:

p.log.restricted <- function(x) {
  
  if (any(x < 0, x > 1)) {
    return(-Inf)
  }
  
  B <- 0.03 # controls 'bananacity'
  -x[1]^2 / 200 - 1 / 2 * (x[2] + B * x[1]^2 - 100 * B)^2
}

res <- MCMCEnsemble(
  p.log.restricted,
  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(res$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.4672 0.2734 0.004991        0.03011
#> b 0.6363 0.2702 0.004933        0.02835
#> 
#> 2. Quantiles for each variable:
#> 
#>      2.5%    25%    50%    75%  97.5%
#> a 0.02576 0.2492 0.4503 0.6889 0.9588
#> b 0.08348 0.4246 0.6867 0.8761 0.9877
plot(res$samples)