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 (or any other non-normal univariate) analogue, so both refuse one. Under "general", every univariate continuous 'lotri' family except wiener() is evaluated, truncated to the parameter's own bounds.

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

# priors need 'lotri' >= 1.0.5
if (utils::packageVersion("lotri") >= "1.0.5") {
  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
#> 
# }