buildMCEM: Builds an MCEM algorithm from a given NIMBLE model

In thirdwing/nimble: Flexible BUGS-compatible system for hierarchical statistical modeling and algorithm development

Description

Takes a nimble model and builds an MCEM algorithm for it. The user must specify which latent nodes are to be integrated out in the E-Step. All other stochastic non-data nodes will be maximized over. If the nodes do not have positive density on the entire real line, then box constraints can be used to enforce this. The M-step is done by a nimble MCMC sampler. The E-step is done by a call to R's optim with method = 'L-BFGS-B'.

Usage

buildMCEM(model, latentNodes, burnIn = 100, mcmcControl = list(adaptInterval
  = 20), boxConstraints = list(), buffer = 10^-6)

Arguments

model	A nimble model
latentNodes	A character vector of the names of the stochastic nodes to integrated out. Names can be expanded, but don't need to be. For example, if the model contains x[1], x[2] and x[3] then one could provide either latentNodes = c('x[1]', 'x[2]', 'x[3]') or latentNodes = 'x'.
burnIn	burn-in used for MCMC sampler in E step
mcmcControl	list passed to MCMCSpec, a nimble function that builds the MCMC sampler. See help(MCMCSpec) for more details
boxConstraints	A list of box constraints for the nodes that will be maximized over. Each constraint is a list in which the first element is a character vector of node names to which the constraint applies and the second element is a vector giving the lower and upper limits. Limits of -Inf or Inf are allowed.
buffer	A buffer amount for extending the boxConstraints. Many functions with boundary constraints will produce NaN or -Inf when parameters are on the boundary. This problem can be prevented by shrinking the boundary a small amount.

Details

buildMCEM calls the NIMBLE compiler to create the MCMC and objective function as nimbleFunctions. If the given model has already been used in compiling other nimbleFunctions, it is possible you will need to create a new copy of the model for buildMCEM to use.

Value

an R function that when called runs the MCEM algorithm. The function returned takes the arguments listed in Runtime Arguments.

See user manual for more

Runtime Arguments

• maxit

  Maximum iterations of the algorithm. Right now, the MCEM runs until maxit, rather than evaluating a more complex stopping criteria

• m1

  number of MCMC samples taken in each iteration before maxit/2

• m2

  number of MCMC samples taken in each iteration after maxit/2

Author(s)

Clifford Anderson-Bergman

Examples

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

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))
pumpModel <- nimbleModel(code = pumpCode, name = 'pump', constants = pumpConsts,
                  data = pumpData, inits = pumpInits)

# Want to maximize alpha and beta (both which must be positive) and integrate over theta
box = list( list(c('alpha','beta'), c(0, Inf)))

pumpMCEM <- buildMCEM(model = pumpModel, latentNodes = 'theta[1:10]',
                       boxConstraints = box)
pumpMCEM(maxit = 40, m1 = 1000, m2 = 5000)

# Could also use latentNodes = 'theta' and buildMCEM would figure out this means 'theta[1:10]'

thirdwing/nimble documentation built on May 31, 2019, 10:41 a.m.