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.12.2 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit https://R-nimble.org.
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)
})
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
Cpump <- compileNimble(pump)
## Compiling
## [Note] This may take a minute.
## [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
box <- list( list(c("alpha","beta"), c(0, Inf))) ## Constraints for the parameters
pumpMCEM <- buildMCEM(model = pump, latentNodes = "theta[1:10]", boxConstraints = box)
## 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.6666402 0.9075246
## Convergence Criterion: 1.001.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 2.
## Current number of MCMC iterations: 1000.
## Parameter Estimates:
## alpha beta
## 0.7794783 1.1415740
## Convergence Criterion: 0.05967816.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 3.
## Current number of MCMC iterations: 1000.
## Parameter Estimates:
## alpha beta
## 0.8120816 1.2321717
## Convergence Criterion: 0.01206093.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## Iteration Number: 4.
## Current number of MCMC iterations: 1000.
## Parameter Estimates:
## alpha beta
## 0.826678 1.280542
## Convergence Criterion: 0.004423875.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## 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: 5.
## Current number of MCMC iterations: 2188.
## Parameter Estimates:
## alpha beta
## 0.8239145 1.2601393
## Convergence Criterion: 0.00149736.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## 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.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## 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.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## 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: 65399.
## Parameter Estimates:
## alpha beta
## 0.823495 1.261259
## Convergence Criterion: 2.188709e-05.
pumpMLE
## alpha beta
## 0.823495 1.261259