This is an old version, view current version.

9.3 1D Integrator

Stan provides a built-in mechanism to perform 1D integration of a function via quadrature methods.

It operates similarly to the algebraic solver and the ordinary differential equations solver in that it allows as an argument a function.

Like both of those utilities, some of the arguments are limited to data only expressions. These expressions must not contain variables other than those declared in the data or transformed data blocks.

9.3.1 Specifying an Integrand as a Function

Performing a 1D integration requires the integrand to be specified somehow. This is done by defining a function in the Stan functions block with the special signature:

real integrand(real x, real xc, real[] theta,
                      real[] x_r, int[] x_i)

The function should return the value of the integrand evaluated at the point x.

The argument of this function are:

• x, the independent variable being integrated over

• xc, a high precision version of the distance from x to the nearest endpoint in a definite integral (for more into see section Precision Loss).

• theta, parameter values used to evaluate the integral

• x_r, data values used to evaluate the integral

• x_i, integer data used to evaluate the integral

Like algebraic solver and the differential equations solver, the 1D integrator separates parameter values, theta, from data values, x_r.

9.3.2 Call to the 1D Integrator

real integrate_1d (function integrand, real a, real b, real[] theta, real[] x_r, int[] x_i)
Integrates the integrand from a to b.

real integrate_1d (function integrand, real a, real b, real[] theta, real[] x_r, int[] x_i, real relative_tolerance)
Integrates the integrand from a to b with the given relative tolerance.

9.3.2.1 Arguments to the 1D Integrator

The arguments to the 1D integrator are as follows:

• integrand: function literal referring to a function specifying the integrand with signature (real, real, real[], real[], int[]):real The arguments represent
  • 1. where integrand is evaluated,
  • 2. distance from evaluation point to integration limit for definite integrals,
  • 3. parameters,
  • 4. real data
  • 5. integer data, and the return value is the integrand evaluated at the given point,
• a: left limit of integration, may be negative infinity, type real,
• b: right limit of integration, may be positive infinity, type real,
• theta: parameters only, type real[],
• x_r: real data only, type real[],
• x_i: integer data only, type int[].

A relative_tolerance argument can optionally be provided for more control over the algorithm:

• relative_tolerance: relative tolerance for the 1d integrator, type real, data only.

9.3.2.2 Return value

The return value for the 1D integrator is a real, the value of the integral.

9.3.2.3 Zero-crossing integrals

For numeric stability, integrals on the (possibly infinite) interval \((a, b)\) that cross zero are split into two integrals, one from \((a, 0)\) and one from \((0, b)\). Each integral is separately integrated to the given relative_tolerance.

9.3.2.4 Precision loss near limits of integration in definite integrals

When integrating certain definite integrals, there can be significant precision loss in evaluating the integrand near the endpoints. This has to do with the breakdown in precision of double precision floating point values when adding or subtracting a small number from a number much larger than it in magnitude (for instance, 1.0 - x). xc (as passed to the integrand) is a high-precision version of the distance between x and the definite integral endpoints and can be used to address this issue. More information (and an example where this is useful) is given in the User’s Guide. For zero crossing integrals, xc will be a high precision version of the distance to the endpoints of the two smaller integrals. For any integral with an endpoint at negative infinity or positive infinity, xc is set to NaN.

9.3.2.5 Algorithmic details

Internally the 1D integrator uses the double-exponential methods in the Boost 1D quadrature library. Boost in turn makes use of quadrature methods developed in (Takahasi and Mori 1974), (Mori 1978), (Bailey, Jeyabalan, and Li 2005), and (Tanaka et al. 2009).

The gradients of the integral are computed in accordance with the Leibniz integral rule. Gradients of the integrand are computed internally with Stan’s automatic differentiation.

References

Takahasi, Hidetosi, and Masatake Mori. 1974. “Double Exponential Formulas for Numerical Integration.” Publications of the Research Institute for Mathematical Sciences 9 (3): 721–41. doi:10.2977/prims/1195192451.

Mori, Masatake. 1978. “An Imt-Type Double Exponential Formula for Numerical Integration.” Publications of the Research Institute for Mathematical Sciences 14 (3): 713–29. doi:10.2977/prims/1195188835.

Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. 2005. “A Comparison of Three High-Precision Quadrature Schemes.” Experiment. Math. 14 (3). A K Peters, Ltd.: 317–29. https://projecteuclid.org:443/euclid.em/1128371757.

Tanaka, Ken’ichiro, Masaaki Sugihara, Kazuo Murota, and Masatake Mori. 2009. “Function Classes for Double Exponential Integration Formulas.” Numerische Mathematik 111 (4): 631–55. doi:10.1007/s00211-008-0195-1.