Numerical integration over (semi-)infinite intervals is easily carried out via changes of variables to map the integrals, e.g., over the finite (0, 1) interval:
integrate(f(x), (x, -inf, +inf)) =
integrate((f((1-t)/t - t/(1-t))) (2t(t-1) + 1) / (t^2(t-1)^2), (t, 0, 1))
integrate(f(x), (x, a, +inf)) = integrate((f(a + (1-t)/t))/t^2, (t, 0, 1))
integrate(f(x), (x, -inf, b)) = integrate((f(b - (1-t)/t))/t^2, (t, 0, 1))
The Jacobian of the variable change is singular at the integration limits that corresponded to an infinte bound, hence this method must be accompanied by integration algorithms that do not evaluate the integrand at the boundaries. For instance, trapezoidal rule would fail to evaluate the integral (in this case a cut-off is required for integral limits, checking convergence to an asymptotic result). Instead, e.g., Gaussian quadrature rules do not necessarily evaluate the integrand at the boundary, so it can be used together with the change of variable above. (Gaussian quadrature actually re-map a finite integral over the interval (-1, 1), but we are not interested here into implementation details of Gaussian rules besides its behavior at the integral limits.)
Let's assume that we have a Racket library (see appendix) providing a procedure (integrate f a b) to integrate f over a finite interval (a, b) with Gaussian quadrature. Then we define procedures to integrate over (semi-)infinte intervals as follows.
;; Integrate over infinite interval (-inf, inf).
(define (integrate-i f)
(define (unit-f t)
(* (f (- (/ (- 1.0 t) t) (/ t (- 1.0 t))))
(/ (+ (* 2.0 t (- t 1.0)) 1.0)
(expt (- t 1.0) 2.0) (* t t))))
(integrate unit-f 0.0 1.0))
;; Integrate with infinite upper limit, (a, +inf).
(define (integrate-iu f a)
(define (unit-f t)
(/ (f (+ a (/ (- 1.0 t) t)))
(* t t)))
(integrate unit-f 0.0 1.0))
;; Integrate with infinite lower limit, (-inf, b).
(define (integrate-il f b)
(define (unit-f t)
(/ (f (- b (/ (- 1.0 t) t)))
(* t t)))
(integrate unit-f 0.0 1.0))
Let's evaluate the integral of a normal distribution with zero mean and unit variance. The expected result is 1 over an infinite range, and 1/2 over semi-infinite ranges (0, +inf) and (-inf, 0).
(define (normal-distribution x)
(/ (exp (- (* 0.5 x x))) (sqrt (* 2.0 pi))))
(integrate-i normal-distribution)
;=> 0.9962262213921168
(integrate-iu normal-distribution 0.0)
;=> 0.49935079589393816
(integrate-il normal-distribution 0.0)
;=> 0.49935079589393816
The change of variable strategy is used, among others, by QUADPACK (see, e.g., its implementation in the GNU Scientific Library). Some library define integrals only within the unit interval, the following change of variables is handy:
integrate(f(x), (x, a, b)) = integrate((f(a + (b-a)t)) (b-a), (t, 0, 1))
It is straightforward to generalize the discussion to the multi-dimensional case.
To run the examples above with Racket the following definitions for the Gauss-Legendre Quadrature rule (taken from Rosetta Code) are required.
#lang racket
;; Gauss-Legendre Quadrature from Rosetta Code.
(define (LegendreP n x)
(let compute ([n n] [Pn-1 x] [Pn-2 1])
(case n
[(0) Pn-2]
[(1) Pn-1]
[else (compute (- n 1)
(/ (- (* (- (* 2 n) 1) x Pn-1)
(* (- n 1) Pn-2)) n)
Pn-1)])))
(define (LegendreP′ n x)
(* (/ n (- (* x x) 1))
(- (* x (LegendreP n x))
(LegendreP (- n 1) x))))
(define (LegendreP-root n i)
; newton-raphson step
(define (newton-step x)
(- x (/ (LegendreP n x) (LegendreP′ n x))))
; initial guess
(define x0 (cos (* pi (/ (- i 1/4) (+ n 1/2)))))
; computation of a root with relative accuracy 1e-15
(if (< (abs x0) 1e-15)
0
(let next ([x′ (newton-step x0)] [x x0])
(if (< (abs (/ (- x′ x) (+ x′ x))) 1e-15)
x′
(next (newton-step x′) x′)))))
(define (Gauss-Legendre-quadrature n)
;; positive roots
(define roots
(for/list ([i (in-range (floor (/ n 2)))])
(LegendreP-root n (+ i 1))))
;; weights for positive roots
(define weights
(for/list ([x (in-list roots)])
(/ 2 (- 1 (sqr x)) (sqr (LegendreP′ n x)))))
;; all roots and weights
(values (append (map - roots)
(if (odd? n) (list 0) '())
(reverse roots))
(append weights
(if (odd? n) (list (/ 2 (sqr (LegendreP′ n 0)))) '())
(reverse weights))))
(define (integrate f a b #:nodes (n 10))
(define m (/ (+ a b) 2))
(define d (/ (- b a) 2))
(define-values [x w] (Gauss-Legendre-quadrature n))
(define (g x) (f (+ m (* d x))))
(* d (+ (apply + (map * w (map g x))))))