aghq {aghq}

Adaptive Gauss-Hermite Quadrature

Description

Normalize the log-posterior distribution using Adaptive Gauss-Hermite Quadrature. This function takes in a function and its gradient and Hessian, and returns a list of information about the normalized posterior, with methods for summarizing and plotting.

Usage

aghq(
  ff,
  k,
  startingvalue,
  transformation = default_transformation(),
  optresults = NULL,
  basegrid = 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.

k

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

startingvalue

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

transformation

Optional. Do the quadrature for parameter theta, but return summaries and plots for parameter g(theta). 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.

basegrid

Optional. Provide an object of class NIGrid from the mvQuad package, representing the base quadrature rule that will be adapted. This is only for users who want more complete control over the quadrature, and is not necessary if you are fine with the default option which basically corresponds to mvQuad::createNIGrid(length(theta),'GHe',k,'product'). Note: the mvQuad functions used within aghq operate on grids in memory, so your basegrid object will be changed after you run aghq.

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.

Details

When k = 1 the AGHQ method is a Laplace approximation, and you should use the aghq::laplace_approximation function, since some of the methods for aghq objects won't work with only one quadrature point. Objects of class laplace have different methods suited to this case. See ?aghq::laplace_approximation.

Value

An object of class aghq which is a list containing elements:

• normalized_posterior: The output of the normalize_logpost function, which itself is a list with elements:

  • nodesandweights: a dataframe containing the nodes and weights for the adaptive quadrature rule, with the un-normalized and normalized log posterior evaluated at the nodes.

  • thegrid: a NIGrid object from the mvQuad package, see ?mvQuad::createNIGrid.

  • lognormconst: the actual result of the quadrature: the log of the normalizing constant of the posterior.

• marginals: a list of the same length as startingvalue of which element j is the result of calling aghq::marginal_posterior with that j. This is a tbl_df/tbl/data.frame containing the normalized log marginal posterior for theta_j evaluated at the original quadrature points. Has columns "thetaj","logpost_normalized","weights", where j is the j you specified.

• optresults: information and results from the optimization of the log posterior, the result of calling aghq::optimize_theta. This a list with elements:

  • ff: the function list that was provided

  • mode: the mode of the log posterior

  • hessian: the hessian of the log posterior at the mode

  • convergence: specific to the optimizer used, a message indicating whether it converged

• control: the control parameters passed

See Also

Other quadrature: get_hessian(), get_log_normconst(), get_mode(), get_nodesandweights(), get_numquadpoints(), get_opt_results(), get_param_dim(), laplace_approximation(), 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))