Estimate log-marginal likelihood using Laplace approximation, but replacing exact calculation of log-determinant with a stochastic approximation
Arguments
- obj
output from
TMB::MakeADFunwhen using penalized likelihood- x
parameter vector used when calculating the Hessian matrix
- k
dimension for Kyrlov subspace
- m
number of probe-vectors to use for approximating average and standard deviation of log-determinant
- seed
if not NULL, then sets the seed. This is helfpul given that the Hutchinson probe vectors are randomly sampled, and comparisons have lower variance using a fixed seed.
Details
This function is only intended when integrating across all parameters, e.g., when supplying a penalized likelihood model with fixed effects mapped off at a prior estimate. For more control over which parameters to estimate, use lanczos_MakeADFun
Examples
# Simulate lognormal-gamma process
set.seed(123)
library(RTMB)
n = 30
n_sum = 3
u = 0 + rnorm(n)
y = rgamma( n, shape = 1/0.5^2, scale = exp(u) * 0.5^2 )
# Fit as GLMM
what = "jnll"
nll = function(p){
nll1 = dnorm(p$u, mean=p$mu, sd=exp(p$logsd), log=TRUE)
nll2 = dgamma(y, shape = 1/exp(2*p$logcv), scale = exp(p$u) * exp(2*p$logcv), log=TRUE)
jnll = -1 * ( sum(nll1) + sum(nll2) )
return(jnll)
}
obj = RTMB::MakeADFun( nll, list(u=u, mu = 0, logsd = 0, logcv = 0), random = "u", silent = TRUE )
opt = nlminb( obj$par, obj$fn, obj$gr )
sdrep = sdreport(obj, bias.correct = TRUE )
H = obj$env$spHess(par = obj$env$last.par.best, random = TRUE)
# Re-do as penalized likelihood
newmap = list(mu = factor(NA), logsd = factor(NA), logcv = factor(NA))
pen = RTMB::MakeADFun( nll, obj$env$parList(), map = newmap, silent = TRUE )
opt_pen = nlminb( pen$par, pen$fn, pen$gr )
# Compare marginal likelihoods
lanczos_nll( pen, k = 10, m = 10 )
#> nll sd_nll
#> 36.46908 0.00000
opt$obj
#> [1] 36.46907