Boost/Math: Numerical integration#

For a general overview, see: Wikipedia [1187], Wikipedia [1194].

Trapezoidal Quadrature#

ctxboost.Trapezoidal(f, a, b, tol, max_refinements)#

where ctx is ctxboost.

Parameters:

f – the function for which the integral is determined.
a – the left border of the integration interval.
b – the right border of the integration interval.
tol – the requested tolerance.
max_refinements – the maximal number of refinements. The default is 12.

Returns:

a tuple (integral, error, L1), where integral is an approximation to the integral, using a trapezoidal quadrature; error is an estimate of the error; and L1 is an estimate of the L1 value.

See also: BoostMath [55], Wikipedia [1212].

If we assume only that the integrand is twice continuously differentiable, we can prove that the error of the composite trapezoidal rule is O(h2). Hence halving the interval only cuts the error by about a fourth, which in turn implies that we must evaluate the function many times before an acceptable accuracy can be achieved.

However, the trapezoidal rule has an astonishing property: If the integrand is periodic, and we integrate it over a period, then the trapezoidal rule converges faster than any power of the step size h. This can be seen by examination of the Euler-Maclaurin summation formula, which relates a definite integral to its trapezoidal sum and error terms proportional to the derivatives of the function at the endpoints and the Bernoulli numbers. If the derivatives at the endpoints are the same or vanish, then the error very nearly vanishes. Hence the trapezoidal rule is essentially optimal for periodic integrands.

Other classes of integrands which are integrated efficiently by this method are the bump functions and bell-shaped integrals over the infinite interval. For details, see Trefethen’s SIAM review.

The routine in Python:

>>> from xlcalcnet import xreal
>>> f = lambda x: 1 / (5 - 4 * xreal.Cos(x))
>>> result2 = xreal.Trapezoidal(f, a, b)
>>> print("res: ", res)

The same routine in Visual Basic:

Function F13 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = 1/(5 - 4*xreal.Cos(x))
    'Console.WriteLine("In F3: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoTrapezoidal()
    Console.WriteLine("Trapezoidal")
    Dim a = 0.0
    Dim b = 2 * xreal.PI
    Dim res1 = xreal.Trapezoidal(AddressOf F13, a, b)
    Console.WriteLine("res1 (integral, error, cond.no.): {0}", res1)
    Console.WriteLine()
End Sub

boost.trapezoidal(f6, a, b, tol)
result:  1.3639170485190606

Gauss-Legendre quadrature#

ctxboost.GaussLegendre(f, a, b)#

where ctx is ctxboost.

Parameters:

f – the function for which the integral is determined.
a – the left border of the integration interval.
b – the right border of the integration interval.

Returns:

a tuple (integral, L1), where integral is an approximation to the integral, using a non-adaptive Gauss-Legendre quadrature; and L1 is an estimate of the L1 value.

See also: Wikipedia [1204], BoostMath [51].

This is intentionally a very simple quadrature routine, it obtains no estimate of the error, and is not adaptive, but is very efficient in simple cases that involve integrating smooth “bell like” functions and functions with rapidly convergent power series.

The Gaussian quadrature routine support both real and complex-valued quadrature.

Internally class gauss has pre-computed tables of abscissa and weights for 7, 15, 20, 25 and 30 points at up to 100-decimal digit precision. That means that using for example, gauss<double, 30>::integrate incurs absolutely zero set-up overhead from computing the abscissa/weight pairs. When using multiprecision types with less than 100 digits of precision, then there is a small initial one time cost, while the abscissa/weight pairs are constructed from strings.

However, for types with higher precision, or numbers of points other than those given above, the abscissa/weight pairs are computed when first needed and then cached for future use, which does incur a noticeable overhead.

The routine in Python:

>>> from xlcalcnet import xreal
>>> a = 0.0; b = 1.0; tol = 1.0E-9
>>> f = lambda x: x * x * xreal.Atan(x)
>>> res = xreal.GaussLegendre(f, a, b)
>>> print("res: ", res)

The same routine in Visual Basic:

Function F14 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = x * x * xreal.Atan(x)
    'Console.WriteLine("In F4: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoGaussLegendre()
    Console.WriteLine("GaussLegendre")
    Dim a = 0.0
    Dim b = 1.0
    Dim res2 = xreal.GaussLegendre(AddressOf F14, a, b)
    Console.WriteLine("res2 (integral, cond.no.): {0}", res2)
    Console.WriteLine()
End Sub

boost.gauss_legendre(f6, a, b, tol)
result:  14.0

Gauss-Kronrod Quadrature#

ctxboost.GaussKronrod(f, a, b, max_depth, tol)#

where ctx is ctxboost.

Parameters:

f – the function for which the integral is determined.
a – the left border of the integration interval.
b – the right border of the integration interval.
max_depth – the maximal number of levels.
tol – the requested tolerance.

Returns:

a tuple (integral, error, L1), where integral is an approximation to the integral, using a Gauss-Kronrod quadrature; error is an estimate of the error; and L1 is an estimate of the L1 value.

See also: Wikipedia [1202], BoostMath [50].

The idea behind Gaussian quadrature is to choose n nodes and weights in such a way that polynomials of order 2n-1 are integrated exactly. However, integration of polynomials is trivial, so it is rarely done via numerical methods. Instead, transcendental and numerically defined functions are integrated via Gaussian quadrature, and the defining problem becomes how to estimate the remainder. Gaussian quadrature alone (without some form of interval splitting) cannot answer this question.

It is possible to compute a Gaussian quadrature of order n and another of order (say) 2n+1, and use the difference as an error estimate. However, this is not optimal, as the zeros of the Legendre polynomials (nodes of the Gaussian quadrature) are never the same for different orders, so 3n+1 function evaluations must be performed. Kronrod considered the problem of how to interleave nodes into a Gaussian quadrature in such a way that all previous function evaluations can be reused, while increasing the order of polynomials that can be integrated exactly. This allows an a posteriori error estimate to be provided while still preserving exponential convergence. Kronrod discovered that by adding n+1 nodes (computed from the zeros of the Legendre-Stieltjes polynomials) to a Gaussian quadrature of order n, he could integrate polynomials of order 3n+1.

The integration routines provided here will perform either adaptive or non-adaptive quadrature, they should be chosen for the integration of smooth functions with no end point singularities.

The number of points specified in the Points template parameter must be an odd number: giving a (N-1)/2 Gauss quadrature as the comparison for error estimation.

Internally class gauss_kronrod has pre-computed tables of abscissa and weights for 15, 31, 41, 51 and 61 Gauss-Kronrod points at up to 100-decimal digit precision. That means that using for example, gauss_kronrod<double, 31>::integrate incurs absolutely zero set-up overhead from computing the abscissa/weight pairs. When using multiprecision types with less than 100 digits of precision, then there is a small initial one time cost, while the abscissa/weight pairs are constructed from strings.

The Gauss-Kronrod quadrature support integrands defined on the real line and returning complex values. In this case, the template argument is the real type, and the complex type is deduced via the return type of the function.

The routine in Python:

>>> from xlcalcnet import xreal
>>> a = 0.0; b = float('+inf')
>>> f = lambda x: xreal.Exp(-x * x / 2)
>>> result3 = xreal.GaussKronrod(f, a, b)
>>> print("result3: ", result3)

The same routine in Visual Basic:

Function F15 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = xreal.Exp(-x*x / 2)
    'Console.WriteLine("In F5: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoGaussKronrod()
    Console.WriteLine("GaussKronrod")
    Dim a = 0.0
    Dim b = xreal.PositiveInfinity
    Dim res3 = xreal.GaussKronrod(AddressOf F15, a, b)
    Console.WriteLine("res3 (integral, error, cond.no.): {0}", res3)
    Console.WriteLine()
End Sub

boost.gauss_kronrod(f6, a, b, max_depth, tol)
result:  14.000000000000002

Double-exponential quadrature: tanh_sinh#

ctxboost.TanhSinh(f, a, b, tol, max_refinements)#

where ctx is ctxboost.

Parameters:

f – the function for which the integral is determined.
a – the left border of the integration interval.
b – the right border of the integration interval.
tol – the requested tolerance.
max_refinements – the maximal number of refinements. The default is 12.

Returns:

a tuple (integral, error, L1, levels), where integral is an approximation to the integral, using a tanh-sinh quadrature; error is an estimate of the error; L1 is an estimate of the L1 value; and levels is the number of levels used.

See also: Wikipedia [1211], BoostMath [48], Okayama et al. [453], Okayama [454].

The tanh-sinh quadrature routine provided by boost is a rapidly convergent numerical integration scheme for holomorphic integrands. By this we mean that the integrand is the restriction to the real line of a complex-differentiable function which is bounded on the interior of the unit disk \(|z| < 1\), so that it lies within the so-called Hardy space. If your integrand obeys these conditions, it can be shown that tanh-sinh integration is optimal, in the sense that it requires the fewest function evaluations for a given accuracy of any quadrature algorithm for a random element from the Hardy space.

Complex integrands are supported.

The routine in Python:

>>> from xlcalcnet import xreal
>>> a = -1.0; b = 1.0
>>> f = lambda x: 5 * x + 7
>>> res = xreal.TanhSinh(f, a, b)
>>> print("res: ", res)

The same routine in Visual Basic:

Function F16 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = 5*x + 7
    'Console.WriteLine("In F6: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoTanhSinh()
    Console.WriteLine("TanhSinh")
    Dim a = -1.0
    Dim b = 1.0
    Dim res4 = xreal.TanhSinh(AddressOf F16, a, b)
    Console.WriteLine("res4 (integral, error, cond.no., level): {0}", res4)
    Console.WriteLine()
End Sub

boost.tanh_sinh(f6, a, b)
result:  14.0

Double-exponential quadrature: sinh_sinh#

ctxboost.SinhSinh(f, tol, max_refinements)#

where ctx is ctxboost.

Parameters:

f – the function for which the integral is determined.
tol – the requested tolerance.
max_refinements – the maximal number of refinements. The default is 9.

Returns:

a tuple (integral, error, L1, levels), where integral is an approximation to the integral, using a sinh-sinh quadrature over the limits of integration \((-\infty, \infty)\); error is an estimate of the error; L1 is an estimate of the L1 value; and levels is the number of levels used.

See also: BoostMath [47], Okayama et al. [453], Okayama [454].

Returns the value of the integral of f over the limits of integration \((-\infty, \infty)\). Complex integrands are supported.

The routine in Python:

>>> from xlcalcnet import xreal
>>> f = lambda x: math.exp(-x * x)
>>> res = xreal.SinhSinh(f)
>>> print("res: ", res)

The same routine in Visual Basic:

Function F17 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = xreal.Exp(-x*x)
    'Console.WriteLine("In F7: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoSinhSinh()
    Console.WriteLine("SinhSinh")
    Dim res5 = xreal.SinhSinh(AddressOf F17)
    Console.WriteLine("res5 (integral, error, cond.no., level): {0}", res5)
    Console.WriteLine()
End Sub

boost.sinh_sinh(f7)
result:  1.7724538509055168

Double-exponential quadrature: exp_sinh#

ctxboost.ExpSinh((f, tol, max_refinements, a, b)#

where ctx is ctxboost.

Parameters:

f – the function for which the integral is determined.
tol – the requested tolerance.
max_refinements – the maximal number of refinements. The default is 12.
a – the left border of the integration interval. The default is 0.
b – the right border of the integration interval. The default is \(+\infty\). Note that a and b cannot both be finite.

Returns:

Returns the value of the integral of f over the limits of integration \((0, \infty)\). Complex integrands are supported.

See also: BoostMath [46], Okayama et al. [453], Okayama [454].

The routine in Python:

>>> from xlcalcnet import xreal
>>> f = lambda x: math.exp(-3 * x)
>>> result3 = xreal.ExpSinh(f)
>>> print("result3: ", result3)

The same routine in Visual Basic:

Function F18 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = xreal.Exp(-3*x)
    'Console.WriteLine("In F8: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoExpSinh()
    Console.WriteLine("ExpSinh")
    Dim res6  = xreal.ExpSinh(AddressOf F18)
    Console.WriteLine("res6 (integral, error, cond.no., level): {0}", res6)
    Console.WriteLine()
End Sub

boost.exp_sinh(f8)
result:  0.33333333333333337

Fourier Integral, Cosine#

ctxboost.OouraCos(f)#

where ctx is ctxboost.

Parameters:: f – the function for which the integral is determined.
Returns:: a tuple (integral, error), where integral is an approximation to the integral, using a Ooura Cos quadrature of the function f in the interval (0, +inf); error is an estimate of the error.

See also: BoostMath [49], Ooura and M. [473], Ooura [474].

The routine in Python:

>>> from xlcalcnet import xreal
>>> f = lambda x: 1 / (x * x + 1)
>>> result3 = xreal.Ooura_Cos(f)
>>> print("result3: ", result3)

The same routine in Visual Basic:

Function F19 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = 1 / (x * x + 1)
    'Console.WriteLine("In F8: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoOoura_Cos()
    Console.WriteLine("Ooura_Cos")
    Dim res7  = xreal.Ooura_Cos(AddressOf F19)
    Console.WriteLine("res7 (integral, error): {0}", res7)
    Console.WriteLine()
End Sub

boost.OouraCos(f3, a, b, get_digits, maxit)
x0:  2.9086877250246185e-09

Fourier Integral, Sine#

ctxboost.OouraSin(f)#

where ctx is ctxboost.

Parameters:: f – the function for which the integral is determined.
Returns:: a tuple (integral, error), where integral is an approximation to the integral, using a Ooura Sin quadrature of the function f in the interval (0, +inf); error is an estimate of the error.

See also: BoostMath [49], Ooura and M. [473], Ooura [474].

The routine in Python:

>>> from xlcalcnet import xreal
>>> f = lambda x: 1 / x
>>> result3 = xreal.Ooura_Sin(f)
>>> print("result3: ", result3)

The same routine in Visual Basic:

Function F20 (ByVal x As xreal) As xreal
    Dim fx As New xreal
    fx = 1 / (x)
    'Console.WriteLine("In F8: x: {0}, f(x): {1}", x, fx)
    Return fx
End Function

Sub DemoOoura_Sin()
    Console.WriteLine("Ooura_Sin")
    Dim res8  = xreal.Ooura_Sin(AddressOf F20)
    Console.WriteLine("res8 (integral, error): {0}", res8)
    Console.WriteLine()
End Sub

boost.OouraCos(f3, a, b, get_digits, maxit)
x0:  2.9086877250246185e-09

Calculating the pdf from the characteristic function#

ctx.pdf_from_cf_boost(x, cf)#

where ctx is dec, mpm, or gmp.

Calculates the pdf as the inverse Fourier transform of its characteristic function.

The PDF of Y is the inverse Fourier transform of its characteristic function,

\[\text{pdf}_X(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{ity} C_X(t) \mathrm{d} y = \frac{1}{\pi} \int_{0}^{\infty} \Re \left ( e^{-itx} C_X(t) \right ) \mathrm{d} t.\]

where \(\Re (z)\) denotes the real part of \(z\). We also have

\[\int_{0}^{\infty} \Re \left ( e^{-itx} C_X(t) \right ) \mathrm{d} t = \int_{0}^{\infty} \Re \left ( C_X(t) \right ) \cos(t x) \mathrm{d} t + \int_{0}^{\infty} \Im \left ( C_X(t) \right ) \sin(t x) \mathrm{d} t.\]

Using the right-hand side of this equation allows for efficient use of the quadrature formula of Fillon, with the half-period \(\omega = \pi/x\).

Calculating the cdf from the characteristic function#

ctx.cdf_from_cf_continuous_boost(x, cf)#

where ctx is dec, mpm, or gmp.

Calculates the cdf from the characteristic function using the procedure of Gil-Pelaez.

Assuming that the characteristic function is absolutely integrable over \((-\infty, \infty)\), Gil-Pelaez derived the following inversion formula which requires integration of a real-valued function only. In particular,

\[\text{cdf}_X(x) = \frac{1}{2} - \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{e^{-itx} C_X(t) - e^{itx} C_X(t)}{it} \mathrm{d} t = \frac{1}{2} - \frac{1}{\pi} \int_{0}^{\infty} \Im \left ( \frac{ e^{-itx} C_X(t)}{t} \right ) \mathrm{d}t.\]

where \(\Im (z)\) denotes the imaginary part of \(z\). We also have

\[\int_{0}^{\infty} \Im \left ( \frac{ e^{-itx} C_X(t)}{t} \right ) \mathrm{d}t = \int_{0}^{\infty} \Im \left ( \frac{C_X(t)}{t} \right ) \cos(t x) \mathrm{d}t - \int_{0}^{\infty} \Re \left ( \frac{C_X(t)}{t} \right ) \sin(t x) \mathrm{d}t.\]

Using the right-hand side of this equation allows for efficient use of the quadrature formula of Fillon, with the half-period \(\omega = \pi/x\).

The python code is currently in Charfun.py

The following code provides a test-suite for the numerical inversion of the characteristic function:

class tests_charfunc(rv_cont):

    def __init__(self, rv2, x = 5, a = 0, b = 2):

        cdf_value = rv2.cdf(x)
        print ("rv2.cdf(x): ", cdf_value)

        rv2.set_x(x)
        plot(rv2.gil_pelaez_imag, [a, b], points=200)
        print

        rv2.set_x(x)
        plot(rv2.gil_pelaez_cos, [a, b], points=200)
        print

        rv2.set_x(x)
        plot(rv2.gil_pelaez_sin, [a, b], points=200)
        print

        rv2.set_x(x)
        I0 = quad(rv2.gil_pelaez_imag, [0, +inf])
        print("Integral: ", I0)
        result0 = 0.5 - I0/pi
        print("result0:", result0 )
        print("diff0:", result0 - cdf_value)


        rv2.set_x(x)
        I1 =quadosc(rv2.gil_pelaez_cos, [0, inf], period=1*pi/x) # half period
        print("I1:", I1 )

        rv2.set_x(x)
        I2 =quadosc(rv2.gil_pelaez_sin, [0, inf], period=1*pi/x) # half period
        print("I2:", I2 )

        I3 = I1 + I2
        print("I3:", I3 )
        print("Int diff:", I3 - I0)
        result3 = 0.5 - I3/pi
        print("result3:", result3 )
        print("diff3:", result3 - cdf_value)

Example: non-central chi-squared distribution

The following code shows the difference between using the generic integration (error: diff0: -6.1218939724378892896e-6) and quadrature for oscillatory functions (error: diff3: 8.4703294725430033907e-22).

mp.dps = 20
print()
print ("Hello mpDistributions local ! ")
print()


a = 0.0
b = 2

n = mpf("5")
x = mpf("10")
rv2 = mpr().chisquare(n)

tests_charfunc(rv2, x, a, b)

This produces the following output (plots are to be added):

Hello mpDistributions local !

rv2.cdf(x):  0.92476475385348782128

Integral:  -1.334418597712864291
result0: 0.92475863195951538339
diff0: -6.1218939724378892896e-6
I1: 0.11817924829245123086
I2: -1.4526170785024453884
I3: -1.3344378302099941575
Int diff: -0.000019232497129866511475
result3: 0.92476475385348782128
diff3: 8.4703294725430033907e-22

Boost/Math: Numerical integration

Contents

Boost/Math: Numerical integration#

Trapezoidal Quadrature#

Gauss-Legendre quadrature#

Gauss-Kronrod Quadrature#

Double-exponential quadrature: tanh_sinh#

Double-exponential quadrature: sinh_sinh#

Double-exponential quadrature: exp_sinh#

Fourier Integral, Cosine#

Fourier Integral, Sine#

Calculating the pdf from the characteristic function#

Calculating the cdf from the characteristic function#