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 overxc, 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 integralx_r, data values used to evaluate the integralx_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[]):realThe arguments represent- where integrand is evaluated,
- distance from evaluation point to integration limit for definite integrals,
- parameters,
- real data
- integer data, and the return value is the integrand evaluated at the given point,
a: left limit of integration, may be negative infinity, typereal,b: right limit of integration, may be positive infinity, typereal,theta: parameters only, typereal[],x_r: real data only, typereal[],x_i: integer data only, typeint[].
A relative_tolerance argument can optionally be provided for more control over the algorithm:
relative_tolerance: relative tolerance for the 1d integrator, typereal, 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.