laplace_check: Assess the accuracy of the Laplace approximation

In RTMBdist: Distributions Compatible with Automatic Differentiation by 'RTMB'

Assess the accuracy of the Laplace approximation

Description

Diagnoses how well the Laplace approximation used by TMB/RTMB approximates the marginal likelihood, via the importance-sampling check of evaluating the integrand relative to the fitted Gaussian.

Usage

laplace_check(obj, nSamples = 1000)

## S3 method for class 'laplace_check'
print(x, digits = 4, ...)

Arguments

obj	A fitted TMB/RTMB object (as returned by MakeADFun) with random effects. The model is assumed to have been optimised, so that obj$env$last.par.best holds the maximum-likelihood parameters and the corresponding mode of the latent variables.
nSamples	Number of Monte Carlo samples drawn from the Gaussian approximation. Default 1000.
x	An object of class "laplace_check".
digits	Number of decimal places for the diagnostic values.
...	Unused.

Details

For fixed parameters \hat\theta, the marginal likelihood is Z = \int f(y\mid x) f(x)\, dx. Writing the importance weight

\tilde\ell(x) = f(y\mid x) f(x) / N(x\mid \hat x, H^{-1}),

the Laplace approximation equals \tilde\ell(\hat x), while E_{x\sim N(\hat x, H^{-1})}[\tilde\ell(x)] = Z exactly. The two coincide if and only if the joint negative log-likelihood is quadratic in the latent variables. The spread of the log-weights therefore measures the departure from this ideal: a standard deviation near zero indicates a near-exact approximation.

Value

An object of class "laplace_check": a list with the log-weights logw, their standard deviation sd_logw, the Laplace and importance-sampling estimates of the log marginal likelihood (lZ_laplace, lZ_is), their difference log_bias, and the relative effective sample size ess_ratio.

Examples

# Chicken weight example; taken from RTMB Introduction vignette
data(ChickWeight)

parameters <- list(
  mua=0,          ## Mean slope
  sda=1,          ## Std of slopes
  mub=0,          ## Mean intercept
  sdb=1,          ## Std of intercepts
  sdeps=1,        ## Residual Std
  a=rep(0, 50),   ## Random slope by chick
  b=rep(0, 50)    ## Random intercept by chick
)

jnll <- function(parms) {
  getAll(ChickWeight, parms, warn=FALSE)
  ## Optional (enables extra RTMB features)
  weight <- OBS(weight)
  ## Initialize joint negative log likelihood
  nll <- 0
  ## Random slopes
  nll <- nll - sum(dnorm(a, mean=mua, sd=sda, log=TRUE))
  ## Random intercepts
  nll <- nll - sum(dnorm(b, mean=mub, sd=sdb, log=TRUE))
  ## Data
  predWeight <- a[Chick] * Time + b[Chick]
  nll <- nll - sum(dnorm(weight, predWeight, sd=sdeps, log=TRUE))
  ## Get predicted weight uncertainties
  ADREPORT(predWeight)
  ## Return
  nll
}

obj <- MakeADFun(jnll, parameters, random=c("a", "b"), silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)

chk <- laplace_check(obj)
chk
# Laplace approximation exact here: linear Gaussian-Gaussian example

RTMBdist documentation built on Oct. 31, 2026, 5:07 p.m.