Short version of my question
Can I only apply the Gauss-Hermite routine with an infinite interval or can I transform the interval?
Long version (reason I'm asking)
I am interested in solving integral equations numerically and one part of this, if I understood it correctly, is the usage of numerical integration routines. So I ventured into this field and implemented the midpoint, the Simpson, and the Gauss-Legendre routine successfully.
Here is my implementation of the Gauss-Legendre routine in R (depending on the package statmod):
gauss_legendre <- function(f, upper, lower, nr_int, ...) {
#1. Compute the nodes; those are the roots of the legendre-polynom
leg_poly <- gauss.quad(nr_int)
x <- leg_poly$nodes
#2. Compute the weights as follows: 2/((1-x^2)Ableitung(P(x))^2)
w <- leg_poly$weights
#3. Compute the integral
x_transformed <- ((x+1)*(upper - lower)/2 + lower)
(upper - lower)/2*sum(w * mapply(f, x_transformed, ...))
}
So basically, the first two parts give me the numbers tabulated here and the third part just computes:
$$\int_a^b f(x)dx = \frac{b-a}{2}\sum\limits_{i=1}^n \omega_i f\left( \frac{b-a}{2}x_i + \frac{a + b}{2} \right)$$
Let's check that it works (note that integrate is just the built-in function in R):
gauss_legendre(function(x) {x^3 + log(x) + sqrt(x)}, 10, 2, 10)
[1] 2528.836
> integrate(function(x) {x^3 + log(x) + sqrt(x)}, 2, 10)
2528.836 with absolute error < 2.8e-11
As you can see, this works with an interval ranging from 2 to 10 by transforming the $x$ values. Given that the interval is well defined, I don't see any problem using that for any interval.
Next, I implemented the Gauss-Hermite routine:
gauss_hermite <- function(f, upper, lower, nr_int, ...) {
#1. Compute the nodes; those are the roots of the Hermite-polynom
her_poly <- gauss.quad(nr_int, "hermite")
x <- her_poly$nodes
#2. Compute the weights according to the function; attention: this already includes the term exp(x^2)
w <- her_poly$weights * exp(x^2)
#3. Compute the integral
sum(w * mapply(f, x))
}
This works for instance when computing the integral of the density function of the standard normal distribution:
gauss_hermite(dnorm, -Inf, Inf, 100)
1
However, I don't get any useful results when the interval doesn't go from $[-\infty; \infty]$. On the other side, in none of the material I read was it mentioned that the Gauss-Hermite approach is limited to that case (I read that it's useful in connection with Normal random variables though). So my question (sorry for the long prolog): can I only apply the Gauss-Hermite routine with an infinite interval or can I transform the interval?
As you probably guessed from the question, I'm not a mathematician, so any help with this would be great.
Best Answer
Yes, there are transformations that transform finite integrals to doubly infinite integrals, like the $\tanh$-substitution
$$\int_a^b f(x) \mathrm dx=\frac{b-a}{2}\int_{-\infty}^\infty f\left(\frac{a+b}{2}+\frac{b-a}{2}\tanh\,u\right)\operatorname{sech}^2 u \,\mathrm du$$
and yes, you can apply Gauss-Hermite quadrature to any suitably transformed doubly infinite integral:
$$\int_{-\infty}^\infty f(v)\mathrm dv\approx\sum_{k=1}^n w_k \exp(x_k^2)f(x_k)$$
BUT
Gauss-Hermite doesn't work very well if the functions you are integrating are not of the form $\exp(-x^2)f(x)$, where $f(x)$ is a function that is well approximated by a polynomial. Remember that the idea of Gaussian quadrature in general is to "factor out" unruly behavior in your integrands, and keep that behavior to the nodes and weights of that quadrature rule. That supposedly leaves you with a remainder that is much more well-behaved for quadrature than the original integrand. In effect, trying to use Gauss-Hermite to integrate things that don't have the factor $\exp(-x^2)$ is like trying to use a toothbrush to scrub a toilet: you'll manage to finish, but it'll take you a while.
To drive the point home, consider the two integrals
$$\begin{align*} \mathcal I_1&=\int_{-\infty}^\infty \operatorname{sech}\,u;\mathrm du=\pi\\ \mathcal I_2&=\int_{-\infty}^\infty \exp(-u^2)\operatorname{sech}\,u\;\mathrm du\approx 1.4790611714495759 \end{align*}$$
Here are a few short Mathematica routines for generating the nodes and weights for Gauss-Hermite quadrature:
Here's a short table of the relative error in using $n$-point Gauss-Hermite quadrature for numerically evaluating these two integrals (done with Mathematica):
$$\begin{array}{c|c|c} n&-\log_{10}\left|\frac{\text{estimate}}{\mathcal I_1}-1\right| &-\log_{10}\left|\frac{\text{estimate}}{\mathcal I_2}-1\right|\\\hline 2&0.58153&1.3066\\ 5&1.0804&2.6953\\ 10&1.6518&4.3750\\ 15&2.0926&5.6971\\ 20&2.4643&6.8235\\ 25&2.7917&7.8217\\ \end{array}$$
Clearly, the convergence of Gauss-Hermite quadrature for $\mathcal I_1$ is rather shabby compared to the relatively quicker convergence for $\mathcal I_2$. Sometimes, you'll get lucky and find a function where Gauss-Hermite performs well even if it does not have an explicit $\exp(-x^2)$ factor, but those things aren't that common.