This is an old version, view current version.

14.2 Integrator Convergence

The integral is performed with the iterative 1D quadrature methods implemented in the Boost library (Agrawal et al. 2017). If the nth estimate of the integral is denoted \(I_n\) and the nth estimate of the norm of the integral is denoted \(|I|_n\), the iteration is terminated when

\[ \frac{{|I_{n + 1} - I_n|}}{{|I|_{n + 1}}} < \text{relative tolerance} \]

The relative_tolerance parameter can be optionally specified as the last argument to integrate_1d and defaults to the square root of the machine epsilon of double precision floating point numbers (that is a relative tolerance of around 1.5e-8).

14.2.1 Zero-crossing Integrals

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)\). This is because the quadrature methods employed internally can have difficulty near zero.

In this case, each integral is separately integrated to the given relative_tolerance.

14.2.2 Avoiding precision loss near limits of integration in definite integrals

If care is not taken, the quadrature can suffer from numerical loss of precision near the endpoints of definite integrals.

For instance, in integrating the pdf of a beta distribution when the values of \(\alpha\) and \(\beta\) are small, most of the probability mass is lumped near zero and one.

The pdf of a beta distribution is proportional to

\[ p(x) \propto x^{\alpha - 1}(1 - x)^{\beta - 1} \]

Normalizing this distribution requires computing the integral of \(p(x)\) from zero to one. In Stan code, the integrand might look like:

real beta(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
  real alpha = theta[1];
  real beta = theta[2];

  return x^(alpha - 1.0) * (1.0 - x)^(beta - 1.0);
}

The issue is that there will be numerical breakdown in the precision of 1.0 - x as x gets close to one. This is because of the limited precision of double precision floating numbers. This integral will fail to converge for values of alpha and beta much less than one.

This is where xc is useful. It is defined, for definite integrals, as a high precision version of the distance from x to the nearest endpoint. To make use of this for the beta integral, the integrand can be re-coded:

real beta(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
  real alpha = theta[1];
  real beta = theta[2];
  real v;

  if(x > 0.5) {
    v = x^(alpha - 1.0) * xc^(beta - 1.0);
  } else {
    v = x^(alpha - 1.0) * (1.0 - x)^(beta - 1.0);
  }

  return v;
}

This version of the integrand will converge for much smaller values of alpha and beta than otherwise possible.

Note, xc is only used for definite integrals. If either the left endpoint is at negative infinity or the right endpoint is at positive infinity, xc will be NaN.

For zero-crossing definite integrals (see section Zero Crossing) the integrals are broken into two pieces (\((a, 0)\) and \((0, b)\) for endpoints \(a < 0\) and \(b > 0\)) and xc is a high precision version of the distance to the limits of each of the two integrals separately. This means xc will be the a high precision version of a - x, x, or b - x, depending on the value of x and the endpoints.

References

Agrawal, Nikhar, Anton Bikineev, Paul Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno Lalande, et al. 2017. “Double-Exponential Quadrature.” https://www.boost.org/doc/libs/1_66_0/libs/math/doc/html/math_toolkit/double_exponential.html.