compute_moment {aghq}

Compute moments

Description

Compute the moment of any function ff using AGHQ.

Usage

compute_moment(obj, ...)

## S3 method for class 'list'
compute_moment(
  obj,
  ff = function(x) 1,
  gg = NULL,
  nn = NULL,
  type = c("raw", "central"),
  method = c("auto", "reuse", "correct"),
  ...
)

## S3 method for class 'aghq'
compute_moment(
  obj,
  ff = function(x) 1,
  gg = NULL,
  nn = NULL,
  type = c("raw", "central"),
  method = c("auto", "reuse", "correct"),
  ...
)

## Default S3 method:
compute_moment(
  obj,
  ff = function(x) 1,
  gg = NULL,
  method = c("auto", "reuse", "correct"),
  ...
)

Arguments

obj

Object of class aghq output by aghq::aghq(). See ?aghq.

...

Used to pass additional argument ff.

ff

Any R function which takes in a numeric vector and returns a numeric vector. Exactly one of ff or gg must be provided. If both are provided, aghq::compute_moment() will use gg, without warning.

gg

The output of, or an object which may be input to aghq::make_moment_function(). See documentation of that function. Exactly one of ff or gg must be provided. If both are provided, aghq::compute_moment() will use gg, without warning.

nn

A numeric scalar. Compute the approximate moment of this order, E(theta^nn|Y). See details. If nn is provided, compute_moment will use it over ff or gg, without warning.

type

Either 'raw' (default) or 'central', see details.

method

Method for computing the quadrature points used to approximate moment. One of 'reuse' (default) or 'correct'. See details. The default SHOULD be 'correct'; it is currently set to 'reuse' to maintain compatibility of results with previous versions. This will be switched in a future major release.

Details

If multiple of nn, gg, and ff are provided, then compute_moment will use nn, gg, or ff, in that order, without warning.

There are several approximations available. The "best" one is obtained by specifying gg and using method = 'correct'. This recomputes the mode and curvature for the function g(theta)posterior(theta), and takes the ratio of the AGHQ approximation to this function to the AGHQ approximation to the marginal likelihood. This obtains the same relative rate of convergence as the AGHQ approximation to the marginal likelihood. It may take a little extra time, and only works for positive, scalar-valued functions g.

method = 'reuse' re-uses the AGHQ adapted points and weights. It's faster than the correct method, because it does not involve any new optimization, it's just a weighted sum. No convergence theory. Seems to work ok in "practice". "Works" for arbitrary g.

Specifying ff instead of gg automatically uses method = 'reuse'. This interface is provided for backwards compatibility mostly. However, one advantage is that it allows for vector-valued functions, in which case it just returns the corresponding vector of approximate moments. Also, it only requires the adapted nodes and weights, not the ability to evaluate the log-posterior and its derivatives, although this is unlikely to be a practical concern.

Specifying a numeric value nn will return the moment E(theta^nn|Y). This automatically does some internal shifting to get the evaluations away from zero, to avoid the inherent problem of multi-modal "posteriors" that occurs when the posterior mode is near zero, and account for the fact that some of the new adapted quadrature points may be negative. So, the actual return value is E(theta^nn + a|Y) - a for a cleverly-chosen value a.

Finally, type='raw' computes raw moments E(g(theta)|Y), where type='central' computes central moments, E(g(theta - E(g(theta)|Y))|Y). See examples.

Value

A numeric vector containing the moment(s) of ff with respect to the joint distribution being approximated using AGHQ.

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)
)
quad <- aghq(funlist2d,7,c(0,0))