Laplace Approximations with adlaplace
This vignette walks through a hierarchical negative-binomial model: simulate data, build an automatic-differentiation (AD) representation as shards, validate derivatives, solve the inner Laplace problem for random effects, and fit outer parameters with a profile likelihood.
Model
Observations follow a negative binomial with mean \(\mu = \exp(\eta)\). We use two equivalent parameterizations: NB size \(r\), and Gamma-mixing standard deviation \(\theta\) (with \(E[Z]=1\), \(\mathrm{sd}(Z)=\theta\) in the Poisson–Gamma construction). They are linked by \(r = 1/\theta^2\).
Size \(r\).
\[ \begin{aligned} f(y;\mu,r) = & \log\Gamma(y + r) - \log\Gamma(r) - \log\Gamma(y + 1) + r \log\left(\frac{r}{r + \mu}\right) + y \log\left(\frac{\mu}{r + \mu}\right) \\ = & \log\Gamma(y + r) - \log\Gamma(r) - \log\Gamma(y + 1) + r \log r - r \log(r + \mu) + y \log(\mu) - y \log(r + \mu) \end{aligned} \]
SD \(\theta\) (\(r = 1/\theta^2\)).
\[ \begin{aligned} f(y;\mu,\theta) = & \log\Gamma\left(y + \theta^{-2}\right) - \log\Gamma\left(\theta^{-2}\right) - \log\Gamma(y + 1) + \theta^{-2} \log\left(\frac{\theta^{-2}}{\theta^{-2} + \mu}\right) + y \log\left(\frac{\mu}{\theta^{-2} + \mu}\right) \\ = & \log\Gamma\left(y + \theta^{-2}\right) - \log\Gamma\left(\theta^{-2}\right) - \log\Gamma(y + 1) - 2 \theta^{-2} \log \theta - \theta^{-2} \log\left(\theta^{-2} + \mu\right) + y \log(\mu) - y \log\left(\theta^{-2} + \mu\right) \end{aligned} \]
The linear predictor is \(\eta = X\beta +
A\gamma\), with two independent random-effect groups in \(A\) and log-standard-deviation parameters
for each group plus a global overdispersion parameter. Densities are
implemented in src/densities.cpp with headers under
src/include/adlaplace/densities.hpp (e.g. negative binomial
observation and hyperparameter terms, plus random_diagonal
for random effects).
Simulated data
We generate 5000 observations, fixed effects, two random-effect
groups, and Poisson–Gamma overdispersion. Outer parameters are stored on
the log-SD scale (transform_theta = TRUE).
Nobs <- if (vignette_full) 5000 else 400
Nrandom1 <- 10
Nrandom2 <- 25
set.seed(0)
X <- Matrix::Matrix(cbind(1, rbinom(Nobs, 1, prob = 0.5)))
Adf <- data.frame(
r1 = sample(Nrandom1, Nobs, replace = TRUE),
r2 = sample(Nrandom2, Nobs, replace = TRUE)
)
AmatList <- list(
r1 = Matrix::sparseMatrix(
i = seq_len(Nobs),
j = Adf$r1
),
r2 = Matrix::sparseMatrix(
i = seq_len(Nobs),
j = Adf$r2
)
)
Amat <- do.call(cbind, AmatList)
beta <- rep(1, ncol(X))
thetaOrig <- c(0.1, 0.1, 0.2) # SD for each random group, then overdispersion
theta <- log(thetaOrig)
gamma <- rnorm(
ncol(Amat),
sd = rep(thetaOrig[1:2], c(Nrandom1, Nrandom2))
)
eta_true <- as.vector(X %*% beta + Amat %*% gamma)
mu_true <- exp(eta_true)
Z <- rgamma(length(mu_true), thetaOrig[3]^(-2), thetaOrig[3]^(-2))
y <- rpois(length(mu_true), mu_true * Z)
config <- list(
beta = beta,
theta = log(thetaOrig),
transform_theta = TRUE,
# Warm-start near the data-generating values (close to the MLE).
gamma = gamma,
obs_groups = adlaplace::obs_groups(
Amat,
num_shards = if (vignette_full) 50L else 8L
),
verbose = FALSE,
package = "adlaplace"
)
n_beta <- length(config$beta)
n_gamma <- length(config$gamma)
n_theta <- length(config$theta)Building AD shards
Each model piece is compiled separately with
ad_pack_ptr(). The full model has four shards: observation
likelihood, hyperparameter (extra) terms, and one shard per
random-effect group. Shards are merged later with
ad_pack().
beta_map and gamma_map map local design
columns of X and A into the global
(beta, gamma, theta) vector: column j of
beta_map has one structural 1 at row
r when local X column j uses
global beta r (and similarly for gamma_map /
A). When ncol(X) > 0,
ncol(beta_map) must equal ncol(X) with exactly
one nonzero per column.
Observations and hyperparameters
The observation shard uses the third \(\theta\) component (global overdispersion). The extra shard contributes the prior on \(\theta\).
data_obs <- adlaplace:::density_data(
y = y,
A = Amat,
X = X,
theta_map = Matrix::sparseMatrix(
i = n_theta,
j = 1L,
dims = c(n_theta, 1L)
),
ad_kind = "observations",
density = "nbinom_obs"
)
data_extra <- adlaplace::density_data(
data_obs,
ad_kind = "parameters",
density = "nbinom_extra"
)
ad_fun_obs <- adlaplace::ad_pack_ptr(data = data_obs, config = config)
ad_fun_extra <- adlaplace::ad_pack_ptr(data = data_extra, config = config)Random effects
Each group gets a diagonal normal prior with its own \(\theta\) index.
model_r1 <- adlaplace:::density_data(
beta_map = n_beta,
gamma_map = Matrix::sparseMatrix(
i = seq_len(Nrandom1),
j = seq_len(Nrandom1),
dims = c(n_gamma, Nrandom1)
),
theta_map = c(1L, n_theta),
ad_kind = "random",
density = "random_diagonal",
precision = rep(1, Nrandom1)
)
ad_fun_r1 <- adlaplace::ad_pack_ptr(data = model_r1, config = config)
model_r2 <- adlaplace:::density_data(
beta_map = n_beta,
gamma_map = Matrix::sparseMatrix(
i = seq.int(Nrandom1 + 1L, length.out = Nrandom2),
j = seq_len(Nrandom2),
dims = c(n_gamma, Nrandom2)
),
theta_map = c(2L, n_theta),
ad_kind = "random",
density = "random_diagonal",
precision = rep(1, Nrandom2)
)
ad_fun_r2 <- adlaplace::ad_pack_ptr(data = model_r2, config = config)Observation likelihood check
Scan the overdispersion parameter while holding \((\beta, \gamma, \theta_{1:2})\) fixed at
the data-generating values. The AD joint log density (observation +
extra shards) should match stats::dnbinom on the full
sample.
xx <- c(beta, gamma, config$theta)
log_theta_true <- config$theta[n_theta]
nb_log_dens <- function(log_theta) {
r <- exp(-2 * log_theta)
sum(stats::dnbinom(y, mu = mu_true, size = r, log = TRUE))
}
ad_nb_log_dens <- function(log_theta) {
xx_scan <- xx
xx_scan[length(xx_scan)] <- log_theta
adlaplace::joint_log_dens(ad_fun_obs, xx_scan, negative = FALSE) +
adlaplace::joint_log_dens(ad_fun_extra, xx_scan, negative = FALSE)
}
seq_log_theta <- log_theta_true + seq(-0.5, 0.5, length.out = 5)
cmp <- data.frame(
theta = exp(seq_log_theta),
stats = vapply(seq_log_theta, nb_log_dens, numeric(1)),
ad = vapply(seq_log_theta, ad_nb_log_dens, numeric(1))
)
cmp$diff <- cmp$ad - cmp$stats
cmp## theta stats ad diff
## 1 0.1213061 -11266.92 -11266.92 -2.122761e-08
## 2 0.1557602 -11254.93 -11254.93 -7.210474e-09
## 3 0.2000000 -11249.51 -11249.51 -1.482113e-08
## 4 0.2568051 -11266.86 -11266.86 -4.954927e-09
## 5 0.3297443 -11335.62 -11335.62 1.293301e-09Combining shards
Shard handles are combined with c(). This uses move
semantics: the individual ad_pack_ptr inputs to
c() are cleared and must not be reused. Keep a plain
combined ptr for serial debug dens/grad/hess, and pass a
clone_ad_pack_ptr() into
ad_pack(..., num_threads > 1) so multi-thread affinity
does not attach to the plain handle.
ad_fun_plain <- c(
ad_fun_obs, ad_fun_extra,
ad_fun_r1, ad_fun_r2
)
ad_pack <- adlaplace::ad_pack(
adlaplace::clone_ad_pack_ptr(ad_fun_plain),
num_threads = vignette_num_threads
)The per-shard log densities sum to the joint log density at a fixed
parameter vector. Serial dens uses ad_fun_plain; parallel
fun_obj_fdfh uses ad_pack.
x_full <- c(config$beta, config$gamma, config$theta)
shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_fun_plain))
by_shard <- vapply(
shards,
function(s) {
adlaplace::joint_log_dens(
ad_fun_plain, x_full,
ad_shards = s, negative = TRUE
)
},
numeric(1)
)
res_fdf <- fun_obj_fdfh(
ad_pack,
c(config$beta, config$theta), config$gamma,
inner = FALSE, verbose = TRUE
)## fun_obj_fdfh: outer, threads = 2, shards = 53, nvars = 40
## ad_pack configured_num_threads = 2
## thread_groups: [0]=27 [1]=26
## env OMP_NUM_THREADS=(unset) OMP_THREAD_LIMIT=(unset) OMP_DYNAMIC=(unset) OMP_MAX_ACTIVE_LEVELS=(unset) OMP_NESTED=(unset) OMP_PROC_BIND=(unset) OMP_PLACES=(unset) OMP_WAIT_POLICY=(unset)
## env MKL_NUM_THREADS=(unset) OPENBLAS_NUM_THREADS=(unset) VECLIB_MAXIMUM_THREADS=(unset)
## omp (serial, pre-CppadParallelScope): in_parallel=0 thread_num=0 num_threads=1 max_threads=4 num_procs=4 dynamic=0 max_active_levels=1
## CppadParallelScope: drain shard eval buffers (entry)...
## CppadParallelScope: setup num_threads=2
## CppadParallelScope: setup done
## outer get_fdfh: before OpenMP threads=2 Nparams=40 hess_size=826 shards=53
## outer get_fdfh parallel_team=2
## outer get_fdfh: after OpenMP api_err_shard=-1
## outer get_fdfh: scatter Hessian...
## outer get_fdfh: scatter done
## CppadParallelScope: drain shard eval buffers (exit)...
## CppadParallelScope: teardown (free_available + free_all)...
## CppadParallelScope: teardown done
## fun_obj_fdfh: done## [1] 9373.413 10786.407 10332.053 7685.215 9148.511## [1] 11217.37 11217.37 11217.37Joint log density derivatives
Gradients and Hessians can be evaluated on all shards or on a subset. Results should agree when the subset is the full model.
log_all <- adlaplace::joint_log_dens(ad_fun_plain, x_full, negative = TRUE)
log_sub <- adlaplace::joint_log_dens(
ad_fun_plain, x_full,
ad_shards = 0:3, negative = TRUE
)
grad_all <- adlaplace::grad(ad_fun_plain, x_full, negative = TRUE)
grad_sub <- adlaplace::grad(
ad_fun_plain, x_full,
ad_shards = 0:3, negative = TRUE
)
hess_all <- adlaplace::hessian(ad_fun_plain, x_full, negative = TRUE)
hess_sub <- adlaplace::hessian(
ad_fun_plain, x_full,
ad_shards = 0:3, negative = TRUE
)
c(
log_diff = log_all - res_fdf$f,
max_abs_grad_diff = max(abs(grad_all - res_fdf$grad)),
max_abs_hess_diff = max(abs(
as.matrix(hess_all) - as.matrix(res_fdf$hessian)
))
)## log_diff max_abs_grad_diff max_abs_hess_diff
## 7.275958e-11 1.164153e-10 2.328306e-10Inner optimization
Given outer parameters \((\beta,
\theta)\), inner_opt() maximizes over \(\gamma\) (equivalently minimizes \(-\log p(y, \gamma \mid \beta, \theta)\)).
The returned solution should track the simulated random effects.
x_outer <- c(config$beta, config$theta)
x_full <- c(config$beta, config$gamma, config$theta)
inner_res <- adlaplace::inner_opt(
parameters = x_outer,
gamma = config$gamma,
ad_pack = ad_pack,
control = list(
maxit = 20L,
report.level = 0,
report.freq = 0,
grad.tol = 1e-4,
cg.tol = 1e-4,
restart.gamma.clamp = Inf
),
deriv = TRUE,
verbose = FALSE
)
adlaplace::joint_log_dens(ad_fun_plain, x_full, negative = TRUE)## [1] 11217.37## [1] 11201.75plot(
inner_res$inner_opt$solution, gamma,
xlab = "estimated", ylab = "true"
)
abline(0, 1, col = "red", lty = 2)
likelihood
log_lik_laplace() runs the inner optimizer and returns
the Laplace approximation to the marginal log likelihood. With
deriv = TRUE, it also returns the profile gradient and
third-order trace terms used in the determinant correction.
OpenMP shard affinity is fixed when you call
ad_pack(..., num_threads = …)
(owner_thread = shard_index %% num_threads);
log_lik_laplace() does not change it. The same assignment
is used for inner optimization, outer derivatives when
deriv = TRUE, and the profile trace
(trace_hinv_t logic, same OpenMP shard groups).
inner_opt() uses one OpenMP/CppAD parallel session for
trust-region inner optimization, outer get_fdfh, and trace
evaluation when deriv = TRUE. joint_log_dens()
/ grad() / hessian() are serial debug APIs:
they error on handles with num_threads > 1. Keep dens on
a plain ad_pack_ptr (or clone_ad_pack_ptr() of
an affined handle, which clears affinity), and use
fun_obj_fdfh / inner_opt /
log_lik_laplace on the multi-thread ad_pack.
For a single-thread ad_pack that also allows dens, use
num_threads = 1L.
For thread-affinity diagnostics, reinstall from a fresh tarball with
PKG_CXXFLAGS="-DDEBUG" R CMD INSTALL adlaplace_*.tar.gz
(use R CMD INSTALL --clean if installing from source).
Configure prints when DEBUG is enabled. With
config$verbose = TRUE, inner_opt prints CppAD
session phase messages (trust region, outer get_fdfh, trace
shards). With deriv = TRUE, trace runs inside that session;
standalone trace_hinv_t() remains for direct use. Each
CppadParallelScope end drains per-shard ADFun taylor /
trace buffers in sequential mode, then flushes per-thread CppAD pools
before the next .Call.
To re-enable CppAD’s thread_alloc cross-thread
return_memory asserts (normally compiled out by R’s
-DNDEBUG), install with
PKG_CXXFLAGS="-DADLAPLACE_CPPAD_ASSERT" (configure appends
-UNDEBUG). See
inst/scripts/check-cppad-threads.sh. Optionally set
ADLAPLACE_HOLD_MEMORY=0 to force
hold_memory(false) as a diagnostic escape hatch (slower;
not a substitute for buffer draining).
res_noderiv <- adlaplace::log_lik_laplace(
x = x_outer,
ad_pack = ad_pack,
config = config,
deriv = FALSE
)## Beginning optimization
##
## iter f nrm_gr status rad CG iter CG result
## 1 11201.748205 3.552441 Continuing 1.000000 9 Reached tolerance
## 2 11201.745371 0.000845 Continuing 1.000000 9 Reached tolerance
##
## Iteration has terminated
## 2 11201.745371 0.000845 Successres_deriv <- adlaplace::log_lik_laplace(
x = x_outer,
ad_pack = ad_pack,
config = config,
deriv = TRUE
)## Beginning optimization
##
## iter f nrm_gr status rad CG iter CG result
## 1 11201.748205 3.552441 Continuing 1.000000 9 Reached tolerance
## 2 11201.745371 0.000845 Continuing 1.000000 9 Reached tolerance
##
## Iteration has terminated
## 2 11201.745371 0.000845 Successdata.frame(
log_lik = c(res_noderiv$log_lik, res_deriv$log_lik),
neg_log_lik = c(res_noderiv$neg_log_lik, res_deriv$neg_log_lik),
inner_fval = c(res_deriv$inner_opt$fval, res_deriv$inner_opt$fval),
row.names = c("deriv = FALSE", "deriv = TRUE")
)## log_lik neg_log_lik inner_fval
## deriv = FALSE -11292.9 11292.9 11201.75
## deriv = TRUE -11292.9 11292.9 11201.75## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -9.3313 -0.8918 -0.3801 -0.5836 -0.3645 6.1659## [1] -48.228900 -87.520712 -2.913501 7.247355 8.198539Outer optimization
outer_fn() and outer_gr() wrap
log_lik_laplace() for use with stats::optim().
A small cache environment stores the current inner \(\gamma\) start.
cache <- new.env(parent = emptyenv())
cache$gamma <- config$gamma
control_inner <- list(
maxit = 30L,
report.level = 0,
report.freq = 0,
restart.gamma.clamp = Inf
)
# Start near the truth (approx. MLE); few outer steps suffice.
x_start <- x_outer + c(0.05, -0.05, 0.1, -0.1, 0.05)
optim_control <- list(
maxit = if (vignette_full) 25L else 5L,
trace = 0, factr = 1e9, pgtol = 1e-3
)
bounds <- list(lower = rep(-5, length(x_outer)), upper = rep(5, length(x_outer)))
common_optim_args <- c(
list(
par = x_start,
fn = adlaplace::outer_fn,
method = "L-BFGS-B",
control = optim_control,
config = config,
ad_pack = ad_pack,
cache = cache,
control_inner = control_inner
),
bounds
)
outer_fit <- do.call(stats::optim, c(common_optim_args, list(gr = adlaplace::outer_gr)))
fit_cmp <- rbind(
optim = outer_fit$par,
start = x_start,
true = c(beta, log(thetaOrig))
)
colnames(fit_cmp) <- c(
paste0("beta[", seq_len(n_beta), "]"),
paste0("log_theta[", seq_len(n_theta), "]")
)
fit_cmp## beta[1] beta[2] log_theta[1] log_theta[2] log_theta[3]
## optim 1.062494 1.011674 -2.242151 -2.580226 -1.634489
## start 1.050000 0.950000 -2.202585 -2.402585 -1.559438
## true 1.000000 1.000000 -2.302585 -2.302585 -1.609438## [1] 0.10622978 0.07575686 0.19505192## [1] 11289.63Derivative checks
Finite-difference checks along one outer coordinate validate the profile gradient, the determinant derivative, and the chain rule through the inner mode \(u(\beta, \theta)\). The finite-difference points should fall on the AD lines.
config$verbose <- FALSE
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 2), mar = c(2, 2, 2, 0), mgp = c(1, 0.5, 0))
x1 <- outer_fit$par
Dpar <- length(x1)
scan_ll <- deriv_scan_laplace(
ad_pack,
x = x1,
index = Dpar,
width = 0.2,
n = 5L,
config = config,
gamma = cache$gamma,
control = control_inner
)
res_mid <- scan_ll$mid_result
plot(scan_ll$grid, scan_ll$value, type = "l",
xlab = expression(theta[3]), ylab = "-log lik")
plot(scan_ll$grid, scan_ll$grad[, Dpar], type = "l",
xlab = expression(theta[3]), ylab = "grad")
points(scan_ll$mid, scan_ll$fd_value, pch = 16)
legend("topright", legend = c("AD", "finite diff"),
lty = c(1, NA), pch = c(NA, 16), bty = "n")
Du <- 2L
plot(scan_ll$grid, scan_ll$u_hat[, Du], type = "l",
xlab = expression(theta[3]), ylab = expression(u[Du]))
plot(scan_ll$grid, scan_ll$dU_index[, Du], type = "l",
xlab = expression(theta[3]),
ylab = expression(d * u[Du] / d * theta[3]))
points(scan_ll$mid, scan_ll$fd_u[, Du], pch = 16)
plot(scan_ll$grid, scan_ll$half_log_det, type = "l",
xlab = expression(theta[3]), ylab = "half log det")
plot(scan_ll$grid, scan_ll$d_det[, Dpar], type = "l",
xlab = expression(theta[3]),
ylab = expression(d / d * theta[3] ~ log ~ det))
points(scan_ll$mid, scan_ll$fd_det, pch = 16)
hessian_outer <- Matrix::forceSymmetric(res_mid$hessian$outer)
hessian_plain <- adlaplace::hessian(
ad_fun_plain,
res_mid$full_parameters,
inner = FALSE,
negative = TRUE
)
max_abs_hess_diff <- max(abs((hessian_outer - hessian_plain)@x), na.rm = TRUE)
max_abs_hess_diff## [1] 2.546585e-10seq_gamma1 <- seq.int(
ad_pack@sizes["beta"] + 1L,
length.out = ad_pack@sizes["gamma"]
)
dU_check <- -solve(
res_mid$hessian$inner,
res_mid$hessian$outer[seq_gamma1, -seq_gamma1]
)
dU_check2 <- -
res_mid$hessian$H_inv %*%
res_mid$hessian$outer[seq_gamma1, -seq_gamma1]
quantile(abs(dU_check - res_mid$gradient$outer$dU))## 0% 25% 50% 75% 100%
## 0.000000e+00 1.666961e-18 5.204170e-18 3.191891e-16 1.554312e-15## 0% 25% 50% 75% 100%
## 0.000000e+00 1.666961e-18 5.204170e-18 3.191891e-16 1.554312e-15## 0% 25% 50% 75% 100%
## 0 0 0 0 0oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 1))
plot(
scan_ll$grid, scan_ll$dU_index[, Du],
type = "l",
xlab = expression(theta[3]),
ylab = expression(d * u[Du] / d * theta[3])
)
points(scan_ll$mid, scan_ll$fd_u[, Du], pch = 16)
points(
rep(res_mid$parameters[Dpar], 2),
c(dU_check[Du, Dpar], dU_check2[Du, Dpar]),
col = c("red", "blue"), cex = c(5, 3)
)
Joint density derivative checks
The same finite-difference idea applies to the full joint log density and its AD gradient and Hessian. Finite-difference points should fall on the AD lines; Hessian panels compare AD Hessian rows to slopes of the AD gradient.
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(4, 2), mar = c(2, 2, 2, 0), mgp = c(1, 0.5, 0))
x_here <- res_mid$full_parameters
Dpar_dens <- length(x_here) - 1L
shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_fun_plain))
scan_j <- deriv_scan_joint(
ad_fun_plain,
x = x_here,
index = Dpar_dens,
width = 0.1,
n = 5L,
ad_shards = shards,
negative = FALSE
)
plot(scan_j$grid, scan_j$grad[, Dpar_dens], type = "l", ylab = "AD gradient")
points(scan_j$mid, scan_j$fd_value, pch = 16)
legend("topright", legend = c("AD", "finite diff"),
lty = c(1, NA), pch = c(NA, 16), bty = "n")
s_par2 <- sort(unique(c(Dpar_dens, 1:3, seq(to = length(x_here), length.out = 5))))
s_par2 <- s_par2[seq(1, min(c(prod(par("mfrow")) - 1L, length(s_par2))))]
for (Dpar2 in s_par2) {
plot(
scan_j$grid, scan_j$hess_row[, Dpar2],
type = "l",
ylab = paste0("H[", Dpar_dens, ",", Dpar2, "]"),
ylim = range(c(scan_j$hess_row[, Dpar2], scan_j$fd_grad[, Dpar2]), na.rm = TRUE)
)
points(scan_j$mid, scan_j$fd_grad[, Dpar2], pch = 16)
}
nicer interface
Code verifying the omnibus interface gives the same results as barebones implementation.
dat <- cbind(
data.frame(y = y),
Adf,
x = X[, 2]
)
md <- model_data(
data = dat,
formula =
nbinom(y, lower = 1e-9, init = 0.2) ~
x + iid(r1, init = 0.1) + iid(r2, init = 0.1)
)
config2 <- list(
transform_theta = TRUE,
obs_groups = obs_groups(
md$term_data$A,
num_shards = 50L
)
)
ad_fun2 <- ad_pack(md, config2, num_threads = vignette_num_threads)
# Serial dens for mean_mvquad: clone clears multi-thread affinity, then
# ad_pack(..., num_threads = 1) rebuilds templates on a dens-safe handle.
ad_fun2_dens <- ad_pack(clone_ad_pack_ptr(ad_fun2@ptr), num_threads = 1L)
md$observations$y@XTp[, 1:5]## 2 x 5 sparse Matrix of class "dgCMatrix"
## 1 2 3 4 5
## x_linear_x 1 . . 1 1
## intercept 1 1 1 1 1## 2 x 5 sparse Matrix of class "dgCMatrix"
##
## [1,] 1 1 1 1 1
## [2,] 1 . . 1 1## term model label init log
## 1 x linear x_linear 0.000000 FALSE
## 2 intercept intercept intercept 0.000000 FALSE
## 3 y nbinom y_nbinom_sd -1.609438 TRUE
## 4 r1 iid r1_iid -2.302585 TRUE
## 5 r2 iid r2_iid -2.302585 TRUE## 3 x 1 sparse Matrix of class "ngCMatrix"
##
## [1,] .
## [2,] .
## [3,] |param_remap <- c(2, 1, 5, 3, 4)
log_lik2 <- log_lik_laplace(
x = res_deriv$parameters[param_remap],
config = config2,
ad_pack = ad_fun2,
deriv = TRUE
)## Beginning optimization
##
## iter f nrm_gr status rad CG iter CG result
## 1 11202.864027 78.480653 Continuing 1.000000 9 Reached tolerance
## 2 11201.745394 0.346539 Continuing 1.000000 9 Reached tolerance
## 3 11201.745371 0.000013 Continuing 1.000000 9 Reached tolerance
##
## Iteration has terminated
## 3 11201.745371 0.000013 Success## [1] -11292.9 -11292.9## [,1] [,2] [,3] [,4] [,5]
## [1,] -87.52071 -48.2289 8.198539 -2.913501 7.247355
## [2,] -87.52071 -48.2289 8.198539 -2.913501 7.247355mean estimation
Posterior mean of the random effects can be approximated with sparse
nested Gauss–Hermite quadrature (mvQuad): nodes are drawn
from the Laplace proposal N(hat(gamma), H^{-1}) and
reweighted by joint_log_dens relative to that Gaussian. Use
a dens-safe clone when the fit handle has
num_threads > 1.
library(mvQuad)
mean_estimate <- mean_mvquad(
parameters = log_lik2$parameters,
mode = log_lik2$inner_opt$solution,
cov = log_lik2$hessian$H_inv,
ad_pack = ad_fun2_dens,
n = 3
)
quantile(mean_estimate - log_lik2$inner_opt$solution)## 0% 25% 50% 75% 100%
## -3.438840e-04 -2.528644e-04 -1.314348e-04 -1.055125e-04 -6.513421e-05