What coordinate transform does scipy.integrate.quad use for infinite bounds?

Question

In SciPy's quad function, it is possible to integrate over an integral with infinite bounds.

import scipy
import numpy as np

def normal_distribution_pdf(x):
    return np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)

result, eps, info = scipy.integrate.quad(normal_distribution_pdf, 0, np.inf, full_output=True)
print(result)

According to the SciPy docs, infinite bounds are handled using a coordinate transformation, using the QUADPACK qagie subroutine.

qagie

handles integration over infinite intervals. The infinite range is mapped onto a finite interval and subsequently the same strategy as in QAGS is applied.

But it is vague about how the infinite range is mapped onto the finite range. Is it something like log(x)? Or another function?

I also examined the QUADPACK documentation. It provides the destination range, [0, 1], but is also vague about how it happens.

I am interested in knowing the specifics of how an infinite range is mapped to the range [0, 1], because I am interested in understanding the limitations of quad().

Answer 1

Integration with one infinite limit

Assuming that A is the lower limit of integration, it uses the following coordinate transformation to transform the range (0, 1] to the infinite range [A, inf).

This results in the integral

(If you'd like to prove that this transformation is correct, it can be done by performing a u-substitution on the right side with X = A + (1 - T) / T to obtain the left side.)

Integration with two infinite limits

This is done identically to the first case, except that for each evaluation of the integrand, the underlying function is evaluated twice: once at F(X) and once at F(-X).

Source: Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2. Pages 64-65.

Implementing in Python

Rather than allowing SciPy to do this coordinate transformation for you, you can do it manually. This has some interesting uses. Specifically, it allows you to look at what the graph of the integrand looks like after this transformation.

Here is a program that plots the integral of the normal PDF before and after this transformation.

import scipy
import numpy as np
import matplotlib.pyplot as plt


def normal_distribution_pdf(x):
    return np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)


def transform_coords_semi_infinite(func, A):
    """Transform function coordinates from (0, 1] to [A, np.inf)"""
    def inner(T):
        return func(A + (1 - T) / T) * (T ** -2)
    return inner


# This program doesn't use this - it's just here for completeness.
def transform_coords_infinite(func):
    """Transform function coordinates from (0, 1] to (-np.inf, np.inf)."""
    def inner(T):
        x = (1 - T) / T
        return (func(x) + func(-x)) * (T ** -2)
    return inner


transformed_pdf = transform_coords_semi_infinite(normal_distribution_pdf, A=0)
# Note: the integral is not defined at 0, because that would correspond to infinity in the transformed coordinates
eps_start = 1e-8
result, eps, info = scipy.integrate.quad(transformed_pdf, eps_start, 1, full_output=True)

xs = np.linspace(0, 5, 100)
plt.title("Normal distribution in X domain")
plt.plot(xs, normal_distribution_pdf(xs))
plt.show()

xs = np.linspace(eps_start, 1, 100)
plt.title("Normal distribution in T domain")
plt.plot(xs, transformed_pdf(xs))

Output:

Another interesting use of this is the points argument to quad(). Normally, this is incompatible with integrating over infinite bounds. However, since the previous example uses the bounds [0, 1], points can be used. However, one must use the coordinate transform t = 1/(1 - a + x) to transform break points from the X domain to the T domain, where a is the lower limit.