laplace_approximation {aghq}

Laplace Approximation

Description

Wrapper function to implement a Laplace approximation to the posterior. A Laplace approximation is AGHQ with k = 1 quadrature points. However, the returned object is of a different class laplace, and a different summary method is given for it. It is especially useful for high-dimensional problems where the curse of dimensionality renders the use of k > 1 quadrature points infeasible. The summary method reflects the fact that the user may be using this for a high-dimensional problem, and no plot method is given, because there isn't anything interesting to plot.

Usage

laplace_approximation(
  ff,
  startingvalue,
  optresults = NULL,
  control = default_control(),
  ...
)

Arguments

ff

A list with three elements:

• fn: function taking argument theta and returning a numeric value representing the log-posterior at theta

• gr: function taking argument theta and returning a numeric vector representing the gradient of the log-posterior at theta

• he: function taking argument theta and returning a numeric matrix representing the hessian of the log-posterior at theta

The user may wish to use numDeriv::grad and/or numDeriv::hessian to obtain these. Alternatively, the user may consider the TMB package. This list is deliberately formatted to match the output of TMB::MakeADFun.

startingvalue

Value to start the optimization. ff$fn(startingvalue), ff$gr(startingvalue), and ff$he(startingvalue) must all return appropriate values without error.

optresults

Optional. A list of the results of the optimization of the log posterior, formatted according to the output of aghq::optimize_theta. The aghq::aghq function handles the optimization for you; passing this list overrides this, and is useful for when you know your optimization is too difficult to be handled by general-purpose software. See the software paper for several examples of this. If you're unsure whether this option is needed for your problem then it probably is not.

control

A list with elements

• method: optimization method to use:

  • 'sparse_trust' (default): trustOptim::trust.optim with method = 'sparse'

  • 'SR1' (default): trustOptim::trust.optim with method = 'SR1'

  • 'trust': trust::trust

  • 'BFGS': optim(...,method = "BFGS")

  Default is 'sparse_trust'.

• optcontrol: optional: a list of control parameters to pass to the internal optimizer you chose. The aghq package uses sensible defaults.

...

Additional arguments to be passed to ff$fn, ff$gr, and ff$he.

Value

An object of class laplace with summary and plot methods. This is simply a list with elements lognormconst containing the log of the approximate normalizing constant, and optresults containing the optimization results formatted the same way as optimize_theta and aghq.

See Also

Other quadrature: aghq(), get_hessian(), get_log_normconst(), get_mode(), get_nodesandweights(), get_numquadpoints(), get_opt_results(), get_param_dim(), marginal_laplace_tmb(), marginal_laplace(), nested_quadrature(), normalize_logpost(), optimize_theta(), plot.aghq(), print.aghqsummary(), print.aghq(), print.laplacesummary(), print.laplace(), print.marginallaplacesummary(), summary.aghq(), summary.laplace(), summary.marginallaplace()

Examples

logfteta2d <- function(eta,y) {
  # eta is now (eta1,eta2)
  # y is now (y1,y2)
  n <- length(y)
  n1 <- ceiling(n/2)
  n2 <- floor(n/2)
  y1 <- y[1:n1]
  y2 <- y[(n1+1):(n1+n2)]
  eta1 <- eta[1]
  eta2 <- eta[2]
  sum(y1) * eta1 - (length(y1) + 1) * exp(eta1) - sum(lgamma(y1+1)) + eta1 +
    sum(y2) * eta2 - (length(y2) + 1) * exp(eta2) - sum(lgamma(y2+1)) + eta2
}
set.seed(84343124)
n1 <- 5
n2 <- 5
n <- n1+n2
y1 <- rpois(n1,5)
y2 <- rpois(n2,5)
objfunc2d <- function(x) logfteta2d(x,c(y1,y2))
funlist2d <- list(
  fn = objfunc2d,
  gr = function(x) numDeriv::grad(objfunc2d,x),
  he = function(x) numDeriv::hessian(objfunc2d,x)
)

thequadrature <- aghq(funlist2d,3,c(0,0))