Nimble Model Construction

library("nimble")
library("lazyeval")

Can we declare the nimbleCode parts in sections? Such as defining the model and priors in separate, reusable code blocks? Here’s a “model” block:

model <- nimbleCode({
for (i in 1:N){
theta[i] ~ dgamma(alpha,beta)
lambda[i] <- theta[i]*t[i]
x[i] ~ dpois(lambda[i])
}
})

We may also wish to define the priors separately, and with specified list of hyperparameters for the priors. Perhaps the need to use quote and to specify each line as a list item, instead of as a block defined with {, is not ideal, but non-standard evaluation is a bit of a wild west still.

hyperparameters <- list(lambda = 1.0, a = 0.1, b = 1.0)

priors <- nimbleCode({
alpha ~ dexp(lambda)
beta ~ dgamma(a,b)
})

Having done so, we can construct a nimbleCode block by passing our chosen hyperparameters to the priors and concatenating the model and priors.

P <- as.expression(sapply(priors, lazyeval::interp, .values = hyperparameters)[-1])
pumpCode <- c(model, P)

The result appears to work as a functional nimbleCode block; though note the resulting expression is not quite identical to the one we would typically write (e.g. without commas):

pumpCode
expression({
for (i in 1:N) {
theta[i] ~ dgamma(alpha, beta)
lambda[i] <- theta[i] * t[i]
x[i] ~ dpois(lambda[i])
}
}, alpha ~ dexp(1), beta ~ dgamma(0.1, 1))

Here we just specify the rest of the model:

pumpConsts <- list(N = 10,
t = c(94.3, 15.7, 62.9, 126, 5.24,
31.4, 1.05, 1.05, 2.1, 10.5))
pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))
pumpInits <- list(alpha = 1,
beta = 1,
theta = rep(0.1, pumpConsts$N)) pump <- nimbleModel(code = pumpCode, name = 'pump', constants = pumpConsts, data = pumpData, inits = pumpInits) and then verify that it runs. Cmodel <- compileNimble(pump) mcmcspec <- configureMCMC(pump, print=FALSE) mcmc <- buildMCMC(mcmcspec) Cmcmc <- compileNimble(mcmc, project = Cmodel) Cmcmc$run(1000)
NULL
samples <- as.data.frame(as.matrix(Cmcmc\$mvSamples))
plot(samples[ , 'alpha'], samples[ , 'beta'], xlab = expression(alpha), ylab = expression(beta))

(The plot shows the same correlation issue that arises without using a block sampler that is illustrated in the manual.)