Running MCEM in NIMBLE

The Monte Carlo expectation maximization (MCEM) algorithm is a method for finding maximum likelihoods in models with latent or missing information. MCEM uses Monte Carlo integration methods to simplify and evaluate potentially high-dimensional integrals.

Here we use the pump model example from BUGS to illustrate how to implement the MCEM algorithm to find maximum likelihood estimates of parameters for statistical models with latent states.

library(nimble, warn.conflicts = FALSE)

## nimble version 0.9.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit http://R-nimble.org.

NIMBLE code for the pump example

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

Set up the model

Assign values to the constants in the model, provide data, initialize the random variables, and then create 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 = 0.1, beta = 0.1,
                  theta = rep(0.1, pumpConsts$N))

## Create the model
pump <- nimbleModel(code = pumpCode, name = "pump", constants = pumpConsts,
                    data = pumpData, inits = pumpInits)

## defining model...

## building model...

## setting data and initial values...

## running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
## checking model sizes and dimensions...
## model building finished.

Compile the model

Cpump <- compileNimble(pump)

## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.

## compilation finished.

Run the MCEM algorithm provided by NIMBLE

box <- list( list(c("alpha","beta"), c(0, Inf))) ## Constraints for the parameters
pumpMCEM <- buildMCEM(model = pump, latentNodes = "theta[1:10]", boxConstraints = box)

## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.

## compilation finished.

## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.

## compilation finished.

## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.

## compilation finished.

## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.

## compilation finished.

pumpMLE <- pumpMCEM$run()

## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 1.
## Current number of MCMC iterations: 1000.
## Parameter Estimates: 
##     alpha      beta 
## 0.6618672 0.9039461 
## Convergence Criterion: 1.001.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 2.
## Current number of MCMC iterations: 1000.
## Parameter Estimates: 
##     alpha      beta 
## 0.7896048 1.1556557 
## Convergence Criterion: 0.07184352.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 3.
## Current number of MCMC iterations: 1000.
## Parameter Estimates: 
##     alpha      beta 
## 0.8079199 1.2370033 
## Convergence Criterion: 0.01101242.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Monte Carlo error too big: increasing MCMC sample size.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 4.
## Current number of MCMC iterations: 1250.
## Parameter Estimates: 
##     alpha      beta 
## 0.8245038 1.2650646 
## Convergence Criterion: 0.002326667.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Monte Carlo error too big: increasing MCMC sample size.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Monte Carlo error too big: increasing MCMC sample size.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 5.
## Current number of MCMC iterations: 2188.
## Parameter Estimates: 
##    alpha     beta 
## 0.819932 1.260650 
## Convergence Criterion: 0.000271179.

pumpMLE

##    alpha     beta 
## 0.819932 1.260650