Maximum likelihood in NIMBLE

We show here that maximum likelihood estimation can be easily done using NIMBLE for models without latent states. This is also possible for models with latent states, but we will need to use numerical integration, the Laplace approximation, or some other techniques (e.g., MCEM).

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.

A simple Gaussian example

First we define a simple Gaussian model and simulate some data for it.

## Model code
gaussianCode <- nimbleCode({
  for (i in 1:N){
    y[i] ~ dnorm(mu, sd = sigma)
  }
  sigma ~ dexp(1)
  mu ~ dnorm(0, sd = 100)
})

## Simulate data
N <- 1000
mu <- 0
sigma <- 5

set.seed(1)
dat <- rnorm(n = N, mean = mu, sd = sigma)

## Setup
gaussianConsts <- list(N = N)
gaussianData <- list(y = dat)
gaussianInits <- list(mu = 0, sigma = 1)

## Create the model
gaussianModel <- nimbleModel(code = gaussianCode, name = "gaussianModel", constants = gaussianConsts,
                             data = gaussianData, inits = gaussianInits)

## 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
CgaussianModel <- compileNimble(gaussianModel)

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

Find the MLEs of parameters mu and sigma

We can take advantage of the R function optim to maximize the log probability of the model with regard to the parameters of interest.

First, we will define a nimbleFunction that sets the values of nodes and returns the log likelihood of those nodes. This function is what optim will be optimizing.

## Calculate the log-likelihood
calcLoglike <- nimbleFunction(
  setup = function(model, wrt, nodes){},
  run = function(parVals = double(1)){
    values(model, wrt) <<- parVals
    ll <- model$calculate(nodes)
    return(ll)
    returnType(double())
  }
)

To obtain maximum likelihood estimates, we simply create an instance of the nimbleFunction and pass it to optim.

## Log-likelihood for the example
wrt <- c("mu", "sigma")
nodes <- gaussianModel$getDependencies(wrt, self = FALSE)
gaussianll <- calcLoglike(gaussianModel, wrt, nodes)
Cgaussianll <- compileNimble(gaussianll, project = gaussianModel)

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

## Maximize the log-likelihood in R
res <- optim(par = c(0, 1), fn = Cgaussianll$run, control = list(fnscale = -1))

## MLEs
res$par

## [1] -0.05838509  5.17082419

## Analytic results for comparison
c(mean(dat), sd(dat) * sqrt((N-1)/N))

## [1] -0.05824071  5.17199126

Our estimate looks good!