| Title: | Fast 'Rcpp' Implementation of Gauss-Hermite Quadrature |
|---|---|
| Description: | Fast, numerically-stable Gauss-Hermite quadrature rules and utility functions for adaptive GH quadrature. See Liu, Q. and Pierce, D. A. (1994) <doi:10.2307/2337136> for a reference on these methods. |
| Authors: | Alexander W Blocker |
| Maintainer: | Alexander W Blocker <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.0.1 |
| Built: | 2024-11-21 06:23:01 UTC |
| Source: | CRAN |
This package provides functions to compute Gauss-Hermite quadrature rules very quickly with a higher degree of numerical stability (tested up to 2000 nodes).
It also provides function for adaptive Gauss-Hermite quadrature, extending Laplace approximations (as in Liu & Pierce 1994).
| Package: | fastGHQuad |
| Type: | Package |
| License: | MIT |
| LazyLoad: | yes |
Alexander W Blocker
Maintainer: Alexander W Blocker <[email protected]>
Golub, G. H. and Welsch, J. H. (1969). Calculation of Gauss Quadrature Rules. Mathematics of Computation 23 (106): 221-230.
Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika, 81(3) 624-629.
gaussHermiteData, aghQuad,
ghQuad
# Get quadrature rule rule <- gaussHermiteData(1000) # Find a normalizing constant g <- function(x) 1/(1+x^2/10)^(11/2) # t distribution with 10 df aghQuad(g, 0, 1.1, rule) # actual is 1/dt(0,10) # Find an expectation g <- function(x) x^2*dt(x,10) # t distribution with 10 df aghQuad(g, 0, 1.1, rule) # actual is 1.25# Get quadrature rule rule <- gaussHermiteData(1000) # Find a normalizing constant g <- function(x) 1/(1+x^2/10)^(11/2) # t distribution with 10 df aghQuad(g, 0, 1.1, rule) # actual is 1/dt(0,10) # Find an expectation g <- function(x) x^2*dt(x,10) # t distribution with 10 df aghQuad(g, 0, 1.1, rule) # actual is 1.25
Convenience function for integration of a scalar function g based upon its Laplace approximation.
aghQuad(g, muHat, sigmaHat, rule, ...)aghQuad(g, muHat, sigmaHat, rule, ...)
g |
Function to integrate with respect to first (scalar) argument |
muHat |
Mode for Laplace approximation |
sigmaHat |
Scale for Laplace approximation ( |
rule |
Gauss-Hermite quadrature rule to use, as produced by
|
... |
Additional arguments for g |
This function approximates
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:
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.
Numeric (scalar) with approximation integral of g from -Inf to Inf.
Alexander W Blocker [email protected]
Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika, 81(3) 624-629.
# 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# 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
Evaluate Hermite polynomial of given degree at given location. This function is provided for demonstration/teaching purposes; this method is not used by gaussHermiteData. It is numerically unstable for high-degree polynomials.
evalHermitePoly(x, n)evalHermitePoly(x, n)
x |
Vector of location(s) at which polynomial will be evaluated |
n |
Degree of Hermite polynomial to compute |
Vector of length(x) values of Hermite polynomial
Alexander W Blocker [email protected]
gaussHermiteData, aghQuad,
ghQuad
Finds real parts of polynomial's roots via eigendecomposition of companion matrix. This method is not used by gaussHermiteData. Only the real parts of each root are retained; this can be useful if the polynomial is known a priori to have all roots real.
findPolyRoots(c)findPolyRoots(c)
c |
Coefficients of polynomial |
Numeric vector containing the real parts of the roots of the polynomial defined by c
Alexander W Blocker [email protected]
gaussHermiteData, aghQuad,
ghQuad
Computes Gauss-Hermite quadrature rule of requested order using Golub-Welsch algorithm. Returns result in list consisting of two entries: x, for nodes, and w, for quadrature weights. This is very fast and numerically stable, using the Golub-Welsch algorithm with specialized eigendecomposition (symmetric tridiagonal) LAPACK routines. It can handle quadrature of order 1000+.
gaussHermiteData(n)gaussHermiteData(n)
n |
Order of Gauss-Hermite rule to compute (number of nodes) |
This function computes the Gauss-Hermite rule of order n using the Golub-Welsch algorithm. All of the actual computation is performed in C/C++ and FORTRAN (via LAPACK). It is numerically-stable and extremely memory-efficient for rules of order 1000+.
A list containing:
x |
the n node positions for the requested rule |
w |
the w quadrature weights for the requested rule |
Alexander W Blocker [email protected]
Golub, G. H. and Welsch, J. H. (1969). Calculation of Gauss Quadrature Rules. Mathematics of Computation 23 (106): 221-230
Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika, 81(3) 624-629.
Convenience function for evaluation of Gauss-Hermite quadrature
ghQuad(f, rule, ...)ghQuad(f, rule, ...)
f |
Function to integrate with respect to first (scalar) argument; this
does not include the weight function |
rule |
Gauss-Hermite quadrature rule to use, as produced by
|
... |
Additional arguments for f |
This function performs classical unidimensional Gauss-Hermite quadrature with the function f using the rule provided; that is, it approximates
by evaluating
Numeric (scalar) with approximation integral of f(x)*exp(-x^2) from -Inf to Inf.
Alexander W Blocker [email protected]
Golub, G. H. and Welsch, J. H. (1969). Calculation of Gauss Quadrature Rules. Mathematics of Computation 23 (106): 221-230.
Liu, Q. and Pierce, D. A. (1994). A Note on Gauss-Hermite Quadrature. Biometrika, 81(3) 624-629.
# Get quadrature rules rule10 <- gaussHermiteData(10) rule100 <- gaussHermiteData(100) # Check that rule is implemented correctly f <- function(x) rep(1,length(x)) if (!isTRUE(all.equal(sqrt(pi), ghQuad(f, rule10), ghQuad(f, rule100)))) { print(ghQuad(f, rule10)) print(ghQuad(f, rule100)) } # These should be 1.772454 f <- function(x) x if (!isTRUE(all.equal(0.0, ghQuad(f, rule10), ghQuad(f, rule100)))) { print(ghQuad(f, rule10)) print(ghQuad(f, rule100)) } # These should be zero# Get quadrature rules rule10 <- gaussHermiteData(10) rule100 <- gaussHermiteData(100) # Check that rule is implemented correctly f <- function(x) rep(1,length(x)) if (!isTRUE(all.equal(sqrt(pi), ghQuad(f, rule10), ghQuad(f, rule100)))) { print(ghQuad(f, rule10)) print(ghQuad(f, rule100)) } # These should be 1.772454 f <- function(x) x if (!isTRUE(all.equal(0.0, ghQuad(f, rule10), ghQuad(f, rule100)))) { print(ghQuad(f, rule10)) print(ghQuad(f, rule100)) } # These should be zero
Calculate coefficients of Hermite polynomial using recursion relation. This function is provided for demonstration/teaching purposes; this method is not used by gaussHermiteData. It is numerically unstable for high-degree polynomials.
hermitePolyCoef(n)hermitePolyCoef(n)
n |
Degree of Hermite polynomial to compute |
Vector of (n+1) coefficients from requested polynomial
Alexander W Blocker [email protected]
gaussHermiteData, aghQuad,
ghQuad