Value and gradient of a model's prior log density, on the natural scale

Every estimation method that implements a prior needs the same thing: turn the prior column on $iniDf (rxode2::rxUiPriors()) into a log-density and its gradient at the current parameter value – as opposed to rxSolve()'s use of the same priors, which draws study-level variability around the model's initial estimate. No Jacobian for an optimizer's own unconstrained-scale reparameterization is applied here; that is specific to the estimation method's own parameterization and is its caller's responsibility to add (chain-ruling gradTheta/gradOmega through d(natural)/d(unconstrained)), the same way adviJacLogDet() (nlmixr2est/src/inner.cpp) already does for the full-Bayes ADVI path.

Usage

rxPriorLogDensity(
  ui,
  theta = NULL,
  omega = NULL,
  method = c("general", "nwpri", "tnpri")
)

Arguments

ui: rxode2 ui model
theta: named numeric vector of current population-parameter values; only the ones a prior is on need to be present
omega: current omega matrix, needed only when a prior touches an omega element or block; NULL otherwise. Must carry every eta in the model (like ui$omega itself, in any name order) even when only one of its blocks has a prior – the underlying C kernel addresses omega positionally (the model's own eta numbering), not by submatrix, so a smaller matrix containing only the referenced block is not enough
method: "general" (default), "nwpri" or "tnpri"; see this file's header comment for what differs. "nwpri" implements NONMEM's own $PRIOR NWPRI omega parameterization; "tnpri" is the same raw-omega om.<eta> treatment as "general" (Monolix's Bayesian-estimation assumption is on the natural parameter, same as here – chol(Omega^-1) is FOCEI's own internal affair, not this kernel's). Neither "nwpri" nor "tnpri" has a Cauchy analogue, so both refuse one.

Value

list with value (scalar log density, summed over every prior term), gradTheta (named numeric, d/dtheta) and gradOmega (a matrix the same dimension as omega, d/dOmega, or NULL when omega was not given). gradOmega is entrywise, treating omega[i, j] and omega[j, i] as independent, the way -Oi %*% Psi %*% Oi-style matrix calculus is usually reported. A caller whose free parameter moves an off-diagonal pair together (a Cholesky or log-Cholesky parameterization, say) needs gradOmega[i, j] + gradOmega[j, i], which is 2 * gradOmega[i, j] for i != j since gradOmega is itself symmetric – but NOT for a diagonal entry, whose own free parameter is gradOmega[i, i] alone; doubling it as well would be wrong

Details

This is a thin R-level convenience shim: the actual math is the pure C++ rxPriorLogDensityEval() (inst/include/rxode2prior.h), exposed through rxode2's C function-pointer table so a downstream package's own C++ objective can call it directly (no R/Rcpp touch, safe from inside an OpenMP-parallel region) instead of round-tripping through R for every evaluation the way this function does.

Author

Matthew L. Fidler

Examples


# \donttest{
one.cmt <- function() {
 ini({
   tka <- 0.45
   tcl <- 1
   tv <- 3.45
   eta.ka ~ 0.6
   add.sd <- c(0, 0.7)
   prior(tka) ~ dnorm(0, 10)
   prior(add.sd) ~ dcauchy(0, 5)
 })
 model({
   ka <- exp(tka + eta.ka)
   cl <- exp(tcl)
   v <- exp(tv)
   d/dt(depot) <- -ka*depot
   d/dt(center) <- ka*depot - cl/v*center
   cp <- center/v
   cp ~ add(add.sd)
 })
}

rxPriorLogDensity(one.cmt, theta=c(tka=0.1, add.sd=0.5))
#>  
#>  
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> $value
#> [1] -5.292545
#> 
#> $gradTheta
#>         tka      add.sd 
#> -0.00100000 -0.03960396 
#> 
#> $gradOmega
#> NULL
#> 
# }