Parallel MCMC Chains

Back to NIMBLE Vignettes

 

The following demonstrates how to execute nchains chains of an MCMC using the parallel package.

The samples are then inspected, and optionally each chain can continue parallel execuation to extend the MCMC run.

 

library(nimble)
library(parallel)

code <- nimbleCode({
    mu ~ dnorm(0, 0.001)
    sigma ~ dunif(0, 1000)
    for(i in 1:N) {
        y[i] ~ dnorm(mu, sd = sigma)
    }
})

set.seed(0)
N <- 10
y <- rnorm(N, 4, 2)

constants <- list(N = N)

data <- list(y = y)

initsFunction <- function() {
    list(mu = rnorm(1),
         sigma = runif(1, 0, 10))
}
    
nchains <- 3   ## number of chains
cl <- makeCluster(nchains, timeout = 86400)    ## 1 day
clusterExport(cl, c('code', 'constants', 'data'))

for(i in seq_along(cl)) {
    set.seed(i)
    inits <- initsFunction()
    clusterExport(cl[i], c('inits', 'i'))
}

samplesList1 <- clusterEvalQ(cl, {
    library(nimble)
    Rmodel <- nimbleModel(code, constants, data, inits)
    conf <- configureMCMC(Rmodel)
    Rmcmc <- buildMCMC(conf)
    Cmodel <- compileNimble(Rmodel)
    Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
    set.seed(i)
    samples <- runMCMC(Cmcmc, 10000)
    return(samples)
})

names(samplesList1) <- paste0('chain', 1:nchains)

 

Inspect the samples in samplesList1:

lapply(samplesList1, dim)

## $chain1
## [1] 10000     2
## 
## $chain2
## [1] 10000     2
## 
## $chain3
## [1] 10000     2

library(basicMCMCplots)
chainsPlot(samplesList1)

 

Now continue running each MCMC chain for an additional 20000 samples:

samplesList2 <- clusterEvalQ(cl, {
    Cmcmc$run(20000, reset = FALSE)
    samples <- as.matrix(Cmcmc$mvSamples)
    return(samples)
})

names(samplesList2) <- paste0('chain', 1:nchains)

lapply(samplesList2, dim)

## $chain1
## [1] 30000     2
## 
## $chain2
## [1] 30000     2
## 
## $chain3
## [1] 30000     2