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 1.1.0 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit https://R-nimble.org.
## 
## Note for advanced users who have written their own MCMC samplers:
##   As of version 0.13.0, NIMBLE's protocol for handling posterior
##   predictive nodes has changed in a way that could affect user-defined
##   samplers in some situations. Please see Section 15.5.1 of the User Manual.

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
##   [Note] Any error reports that follow may simply reflect missing values in model variables.

## Checking model sizes and dimensions

Compile the model

Cpump <- compileNimble(pump)

## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

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)

##   [Warning] No samplers assigned for 2 nodes, use conf$getUnsampledNodes() for node names.

## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

pumpMLE <- pumpMCEM$run()

## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 1.
## Current number of MCMC iterations: 1000.
## Parameter Estimates: 
##     alpha      beta 
## 0.6654571 0.9045308 
## Convergence Criterion: 1.001.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 2.
## Current number of MCMC iterations: 1000.
## Parameter Estimates: 
##     alpha      beta 
## 0.7779274 1.1255696 
## Convergence Criterion: 0.05836946.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 3.
## Current number of MCMC iterations: 1000.
## Parameter Estimates: 
##    alpha     beta 
## 0.823343 1.243904 
## Convergence Criterion: 0.015189.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 4.
## Current number of MCMC iterations: 1000.
## Parameter Estimates: 
##     alpha      beta 
## 0.8280817 1.2785976 
## Convergence Criterion: 0.003229742.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

## 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: 1625.
## Parameter Estimates: 
##    alpha     beta 
## 0.818744 1.255177 
## Convergence Criterion: 0.001503264.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

## Monte Carlo error too big: increasing MCMC sample size.

## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

## Monte Carlo error too big: increasing MCMC sample size.

## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

## Monte Carlo error too big: increasing MCMC sample size.

## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 6.
## Current number of MCMC iterations: 4298.
## Parameter Estimates: 
##     alpha      beta 
## 0.8227655 1.2644472 
## Convergence Criterion: 0.0002503726.

pumpMLE

##     alpha      beta 
## 0.8227655 1.2644472