marginal_laplace {aghq}

Marginal Laplace approximation

Description

Implement the marginal Laplace approximation of Tierney and Kadane (1986) for finding the marginal posterior (theta | Y) from an unnormalized joint posterior (W,theta,Y) where W is high dimensional and theta is low dimensional. See the AGHQ software paper for a detailed example, or Stringer et. al. (2020).

Usage

marginal_laplace(
  ff,
  k,
  startingvalue,
  transformation = default_transformation(),
  optresults = NULL,
  control = default_control_marglaplace(),
  ...
)

Arguments

ff

A function list similar to that required by aghq. However, each function now takes arguments W and theta. Explicitly, this is a list containing elements:

• fn: function taking arguments W and theta and returning a numeric value representing the log-joint posterior at W,theta

• gr: function taking arguments W and theta and returning a numeric vector representing the gradient with respect to W of the log-joint posterior at W,theta

• he: function taking arguments W and theta and returning a numeric matrix representing the hessian with respect to W of the log-joint posterior at W,theta

k

Integer, the number of quadrature points to use. I suggest at least 3. k = 1 corresponds to a Laplace approximation.

startingvalue

A list with elements W and theta, which are numeric vectors to start the optimizations for each variable. If you're using this method then the log-joint posterior should be concave and these optimizations should not be sensitive to starting values.

transformation

Optional. Do the quadrature for parameter theta, but return summaries and plots for parameter g(theta). This applies to the theta parameters only, not the W parameters. transformation is either: a) an aghqtrans object returned by aghq::make_transformation, or b) a list that will be passed to that function internally. See ?aghq::make_transformation for details.

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 for the theta optimization:

  • 'sparse_trust' (default): trustOptim::trust.optim

  • 'sparse': trust::trust

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

• inner_method: optimization method to use for the W optimization; same options as for method

Default inner_method is 'sparse_trust' and default method is 'BFGS'.

...

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

Value

If k > 1, an object of class marginallaplace, which includes the result of calling aghq::aghq on the Laplace approximation of (theta|Y), .... See software paper for full details. If k = 1, an object of class laplace which is the result of calling aghq::laplace_approximation on the Laplace approximation of (theta|Y).

See Also

Other quadrature: aghq(), get_hessian(), get_log_normconst(), get_mode(), get_nodesandweights(), get_numquadpoints(), get_opt_results(), get_param_dim(), laplace_approximation(), marginal_laplace_tmb(), 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))
objfunc2dmarg <- function(W,theta) objfunc2d(c(W,theta))
objfunc2dmarggr <- function(W,theta) {
  fn <- function(W) objfunc2dmarg(W,theta)
  numDeriv::grad(fn,W)
}
objfunc2dmarghe <- function(W,theta) {
  fn <- function(W) objfunc2dmarg(W,theta)
  numDeriv::hessian(fn,W)
}

funlist2dmarg <- list(
  fn = objfunc2dmarg,
  gr = objfunc2dmarggr,
  he = objfunc2dmarghe
)