[R] Gaussian Adaptive Quadrature

[Earl F. Glynn]

"Ravi Varadhan" <rvaradhan at jhmi.edu> wrote in message 
news:000001c76bcd$5691c550$7c94100a at win.ad.jhu.edu...
> The terminology "Gaussian quadrature" is not restricted to Gauss-Hermite
> quadrature (where exp(-x^2) is the weight function), but applies more
> broadly to Gauss-Legendre, Gauss-Laguerre, etc., where the abscissa are
> chosen from Legendre, Laguerre polynomials.

Also, the functional form of the integrand and limits of integration vary a 
bit with the different Gaussian Quadrature methods.  With Gaussian 
Quadrature the integral from a to b of f(x) dx is approximated as the sum of 
certain weights multiplied by the function evaluated at certain points 
related to the roots of the orthogonal polynomials.

Gauss-Legendre is perhaps the most general, since it works well with most 
functions over a fixed, finite, interval, normally [-1,1].  The Legendre 
polynomials are orthogonal on the interval [-1,1] with respect to a 
weighting function w(x) = 1.  A simple transformation allows the integration 
interval to be any finite interval, [a,b].

Gauss-Laguerre assumes a weighting function w(x) = exp(-x) in the integrand, 
and an interval of integration from 0 to infinity.

Gauss-Hermite assumes a weighting function w(x) = exp(-x^2) in the 
integrand, and an interval of integration from -infinity to infinity.

Applied Numerical Methods by Carnahan et al provides good details and 
examples (but in FORTRAN).


One "adaptive" approach that can be used with Gaussian Quadrature is to use 
a different number of  terms to evaluate the integral.  To save computation 
time, you can use fewer terms.  This file gives the weights needed for 
various N-point Gauss-Legendre quadrature approximation:
http://www.math.ntnu.no/num/nnm/Program/Numlibc/gauss_co.c

Some years ago on a project we found that 2-point Gaussian Quadrature gave 
us an answer that was "good enough," and obviously was quite fast with so 
few function evaluations.  For your problem you might try 2-point to 
15-point quadrature to see if you get the desired accuracy.   I've always 
used pre-computed polynomial roots and weights for the various N-point 
formulas.  I'm not sure how gauss.quad in library(statmod) gets these 
values. It wasn't obvious to me from a quick look at the source code.

Another  "adaptive" Gaussian approach might break a single integral up into 
a number of other integrals.  One could even use different N-point formulas 
over different intervals, using lower N for "smoother" areas, and larger N 
if a function wasn't so well-behaved.

Some other good links:

Gauss-Legendre Quadrature
http://math.fullerton.edu/mathews/n2003/gaussianquad/GaussianQuadBib/Links/GaussianQuadBib_lnk_1.html
http://mathworld.wolfram.com/Legendre-GaussQuadrature.html

efg

Earl F. Glynn
Scientific Programmer
Stowers Institute for Medical Research