marginal_laplace: Marginal Laplace approximation
In awstringer1/aghq: Adaptive Gauss Hermite Quadrature for Bayesian Inference

marginal_laplace

R Documentation

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

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
)

awstringer1/aghq documentation built on March 4, 2024, 4:48 p.m.