nimble (version 0.6-6)

buildMCEM: Builds an MCEM algorithm from a given NIMBLE model

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' if the nodes are constrainted, or method = 'BFGS' if the nodes are unconstrained.

Usage

buildMCEM(model, latentNodes, burnIn = 500, mcmcControl = list(adaptInterval
  = 100), boxConstraints = list(), buffer = 10^-6, alpha = 0.25,
  beta = 0.25, gamma = 0.05, C = 0.001, numReps = 300, verbose = TRUE)

Value

an R list with two elements:

• run A function that when called runs the MCEM algorithm. This function takes the arguments listed in run Arguments below.

• estimateCov An EXPERIMENTAL function that when called estimates the asymptotic covariance of the parameters. The covariance is estimated using the method of Louis (1982). This function takes the arguments listed in estimateCov Arguments below.

<code>run</code> Arguments

• initM starting number of iterations for the algorithm.

<code>estimateCov</code> Arguments

• MLEs named vector of MLE values. Must have a named MLE value for each stochastic, non-data, non-latent node. If the run() method has alread been called, MLEs do not need to be provided.

• useExistingSamples logical argument. If TRUE and the run() method has previously been called, the covariance estimation will use MCMC samples from the last step of the MCEM algorithm. Otherwise, an MCMC algorithm will be run for 10,000 iterations, and those samples will be used. Defaults to FALSE.

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. Uses an ascent-based MCEM algorithm, which includes rules for automatically increasing the number of MC samples as iterations increase, and for determining when convergence has been reached. Constraints for parameter values can be provided. If contstraints are not provided, they will be automatically determined by NIMBLE.

References

Caffo, Brian S., Wolfgang Jank, and Galin L. Jones (2005). Ascent-based Monte Carlo expectation-maximization. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2), 235-251.

Louis, Thomas A (1982). Finding the Observed Information Matrix When Using the EM Algorithm. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 44(2), 226-233.

Examples

# NOT RUN {
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)
MLEs <- pumpMCEM$run(initM = 1000)
cov <- pumpMCEM$estimateCov()
# }
# NOT RUN {
# Could also use latentNodes = 'theta' and buildMCEM() would figure out this means 'theta[1:10]'

# }