aghQuad() R function from [fastGHQuad]

Adaptive Gauss-Hermite quadrature using Laplace approximation

Convenience function for integration of a scalar function g based upon its Laplace approximation.


aghQuad(g, muHat, sigmaHat, rule, ...)

Arguments

g: Function to integrate with respect to first (scalar) argument
muHat: Mode for Laplace approximation
sigmaHat: Scale for Laplace approximation (sqrt(-1/H), where H is the second derivative of g at muHat)
rule: Gauss-Hermite quadrature rule to use, as produced by gaussHermiteData
...: Additional arguments for g

Returns

Numeric (scalar) with approximation integral of g from -Inf to Inf.

Details

This function approximates

\int_{-\infty}^{\infty} g(x) \, dxintegral( g(x), -Inf, Inf)

using the method of Liu & Pierce (1994). This technique uses a Gaussian approximation of g (or the distribution component of g, if an expectation is desired) to "focus" quadrature around the high-density region of the distribution. Formally, it evaluates:

\sqrt{2} \hat{\sigma} \sum_i w_i \exp(x_i^2) g(\hat{\mu} + \sqrt{2}sqrt(2) * sigmaHat * sum( w * exp(x^2) * g(muHat + sqrt(2) * sigmaHat * x))

\hat{\sigma} x_i)sqrt(2) * sigmaHat * sum( w * exp(x^2) * g(muHat+ sqrt(2) * sigmaHat * x))

where x and w come from the given rule.

This method can, in many cases (where the Gaussian approximation is reasonably good), achieve better results with 10-100 quadrature points than with 1e6 or more draws for Monte Carlo integration. It is particularly useful for obtaining marginal likelihoods (or posteriors) in hierarchical and multilevel models --- where conditional independence allows for unidimensional integration, adaptive Gauss-Hermite quadrature is often extremely effective.

Examples


# Get quadrature rules
rule10  <- gaussHermiteData(10)
rule100 <- gaussHermiteData(100)

# Estimating normalizing constants
g <- function(x) 1/(1+x^2/10)^(11/2) # t distribution with 10 df
aghQuad(g, 0, 1.1,  rule10)
aghQuad(g, 0, 1.1,  rule100)
# actual is
1/dt(0,10)

# Can work well even when the approximation is not exact
g <- function(x) exp(-abs(x)) # Laplace distribution
aghQuad(g, 0, 2,  rule10)
aghQuad(g, 0, 2,  rule100)
# actual is 2

# Estimating expectations
# Variances for the previous two distributions
g <- function(x) x^2*dt(x,10) # t distribution with 10 df
aghQuad(g, 0, 1.1,  rule10)
aghQuad(g, 0, 1.1,  rule100)
# actual is 1.25

# Can work well even when the approximation is not exact
g <- function(x) x^2*exp(-abs(x))/2 # Laplace distribution
aghQuad(g, 0, 2,  rule10)
aghQuad(g, 0, 2,  rule100)
# actual is 2

Author(s)

Alexander W Blocker ablocker@gmail.com

References

Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika, 81(3) 624-629.

aghQuad function