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.

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)

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.121124e-08
## 2 0.1557602 -11254.93 -11254.93 -7.215931e-09
## 3 0.2000000 -11249.51 -11249.51 -1.481203e-08
## 4 0.2568051 -11266.86 -11266.86 -5.011316e-09
## 5 0.3297443 -11335.62 -11335.62  1.267836e-09

Combining 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=10 num_procs=10 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

by_shard[seq(1, min(c(5, length(by_shard))))]

## [1]  9373.413 10786.407 10332.053  7685.215  9148.511

c(
  sum(by_shard),
  adlaplace::joint_log_dens(ad_fun_plain, x_full, negative = TRUE),
  res_fdf$f
)

## [1] 11217.37 11217.37 11217.37

Joint 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      9.313226e-10

Inner 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

inner_res$inner_opt$fval

## [1] 11201.75

plot(
  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                    Success

res_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                    Success

data.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

summary(res_deriv$gradient$outer$trace3)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -9.3313 -0.8918 -0.3801 -0.5836 -0.3645  6.1659

res_deriv$deriv$d_neg_log_lik

## [1] -48.228900 -87.520712  -2.913501   7.247355   8.198539

Outer 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

exp(outer_fit$par[n_beta + seq_len(n_theta)])

## [1] 0.10622978 0.07575686 0.19505192

outer_fit$value

## [1] 11289.63

Derivative 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] 4.656613e-10

par(oldpar)

seq_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 2.168404e-18 6.288373e-18 3.122502e-16 1.859624e-15

quantile(abs(dU_check - dU_check2))

##           0%          25%          50%          75%         100% 
## 0.000000e+00 2.168404e-18 6.288373e-18 3.122502e-16 1.859624e-15

quantile(abs(dU_check2 - res_mid$gradient$outer$dU))

##   0%  25%  50%  75% 100% 
##    0    0    0    0    0

oldpar <- 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)
)

par(oldpar)

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)
}

par(oldpar)

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

data_obs@XTp[, 1:5]

## 2 x 5 sparse Matrix of class "dgCMatrix"
##               
## [1,] 1 1 1 1 1
## [2,] 1 . . 1 1

md$term_data$info$parameters[, c("term", "model", "label", "init", "log")]

##        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

data_obs@theta_map

## 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

c(log_lik2$log_lik, res_deriv$log_lik)

## [1] -11292.9 -11292.9

rbind(log_lik2$deriv$d_neg_log_lik, res_deriv$deriv$d_neg_log_lik[param_remap])

##           [,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.247355

mean 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.513422e-05

plot(log_lik2$inner_opt$solution, mean_estimate)
abline(0, 1, col = "red", lty = 2)