add integrate_1d_gauss_kronrod

[avehtari]

Closes #3000.

Summary

Adds integrate_1d_gauss_kronrod, a 1-D quadrature function that wraps
Boost's adaptive Gauss-Kronrod (boost/math/quadrature/gauss_kronrod.hpp,
already in the vendored Boost 1.87.0) in the (almost) same user-facing API shape
as integrate_1d. The two functions cover disjoint weaknesses and live
side by side:

• integrate_1d (tanh_sinh / exp_sinh / sinh_sinh) excels at integrands
  with algebraic or logarithmic endpoint singularities; its
  double-exponential nodes concentrate near the endpoints.
• integrate_1d_gauss_kronrod (K21 with adaptive bisection) excels at
  smooth integrands and at integrands with peaks off-centre in the
  integration interval — exactly the case where double-exponential
  quadrature undersamples and integrate_1d can return silently-zero
  or wildly-imprecise values.

The API is the same as with integrate_1d except adding abs_tol and max_depth arguments.

My use case is integrals of functions which are product of normal
density and a smooth function. Due to the light tails of the normal
density this functions goes to zero in tails and Gauss-Kronrod excels
compare to double exponential. I really needed this for some
experiments I'm running, and abs_tol was also required to get rid of
numerical problems.

Public API

C++ (prim/rev/fwd autodiff layers):

double integrate_1d_gauss_kronrod(
    const F& f, double a, double b,
    const std::vector<double>& theta,
    const std::vector<double>& x_r,
    const std::vector<int>& x_i,
    std::ostream* msgs,
    double relative_tolerance = std::sqrt(EPSILON),
    double absolute_tolerance = 0.0,
    int max_depth = 15);

Stan language — four overloads, parameter additions trailing:

integrate_1d_gauss_kronrod(f, a, b, theta, x_r, x_i);
integrate_1d_gauss_kronrod(f, a, b, theta, x_r, x_i, rel_tol);
integrate_1d_gauss_kronrod(f, a, b, theta, x_r, x_i, rel_tol, abs_tol);
integrate_1d_gauss_kronrod(f, a, b, theta, x_r, x_i, rel_tol, abs_tol, max_depth);

Convergence test is QUADPACK-style mixed:

error <= max(relative_tolerance * L1, absolute_tolerance)

With absolute_tolerance = 0 (default) this reduces to the strict
pure-relative behaviour of integrate_1d.

User integrand functor signature is unchanged from integrate_1d
((x, xc, theta, x_r, x_i, msgs) -> real); xc is unused under
Gauss-Kronrod (always passed as NaN).

Design notes

• K21 default Kronrod order. Matches scipy.integrate.quad,
  MATLAB's quad, QUADPACK's dqag with default key. Polynomial
  exactness up to degree 31 covers the Gaussian-times-likelihood
  shapes that dominate Stan use. Higher N tightens node spacing for
  hard peaks but at 1.5x-3x cost on smooth integrands; we chose to
  address robustness on hard cases through absolute_tolerance and
  informed bound choice instead.
• absolute_tolerance argument. Boost's gauss_kronrod accumulates
  error += error_local across adaptive subdivisions, with each
  error_local carrying a floor 2 * eps * |K_local|. The
  accumulated floor scales as ~2^max_depth * eps * L1, about
  7e-12 * L1 at the default max_depth = 15. For integrals whose
  L1 falls below ~1e-8 (e.g. nested calls where the outer
  evaluates the integrand deep in the tail of a Gaussian factor),
  the strict pure-relative test measures noise against noise and
  throws spuriously. QUADPACK-style mixed convergence is the
  established fix.
• max_depth exposed. Default 15 matches Boost; user can raise it
  for hard adaptive cases without re-implementing the wrapper.
• Deliberate weakness on endpoint singularities. Documented in the
  test suite (a endpoint_singularity_throws test asserts that GK
  throws on 1/sqrt(x), 1/sqrt(1-x), and beta-with-small-shapes
  integrands). Users with these integrands should keep using
  integrate_1d.

Implementation layout

stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp     prim worker + variadic + legacy wrapper
stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp      rev specialisation, nested-autodiff gradients
stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp      fwd specialisation via finite_diff
stan/math/{prim,rev,fwd}/functor.hpp                      add the new include line
test/unit/math/{prim,rev,mix}/functor/                    25 tests, all passing

The user-facing functor signature is adapted through the existing
integrate_1d_adapter (no new adapter needed). stanc3 changes will be in a
companion PR.

Testing

• 25 unit tests across prim (7), rev (17), and mix (1) all passing.
  Tests are ported from integrate_1d's counterparts with adaptations
  to drop the xc-based functors; new tests cover the
  absolute_tolerance and max_depth arguments and document the
  endpoint-singularity weakness.
• A small cmdstan end-to-end model exercises all four Stan-language
  overloads on integrals with known closed forms; rev autodiff is
  verified by integrating a Normal PDF with a parameters-block mu
  across 200 MCMC draws (returns exactly 1 each time).
• A nested-integration stress test (Student-t LMM marginal likelihood
  over a 2-D random-effects vector, 144 observations, 4 chains x 100
  draws) produces ELPD estimates that agree with independent bridge
  sampling to 0.05 nats per fold at the 99th percentile. For the
  same case, the existing integrate_1d was silently returning zeros.

Possible future work

The Kronrod order N is currently fixed at 21 in the wrapper. Boost
provides compile-time tables for N in {15, 21, 31, 41, 51, 61}, and
exposing N as an optional data int parameter via a small
runtime-to-compile-time dispatch (a switch(N) over the six valid
values) is a clean extension that doesn't break the current overload
set. Left out from this PR

Companion PRs

• stanc3: four new integrate_1d_gauss_kronrod signatures in
  Stan_math_signatures.ml and an |-alternative pattern extension
  in Lower_expr.ml's msgs-injection clause.
  (deferred until this PR is approved).
• stan-dev/docs: function reference entry alongside integrate_1d
  (deferred until this PR is approved).

Checklist

• Copyright holder: (fill in copyright holder information)

Stan Development Team. The code is duplication of existing integrate_1d code and there is no additional creativity.

The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

AI Use Disclosure

The duplication of the existing integrate_1d code and tests was made using Claude AI Agent. The actual quadrature algorithm is in Boost library, and the code is just a wrapper to call Boost.

[WardBrian]

Would it make sense to add a suffixed alias of the existing integrate_1d as part of this? I think we have done similar in the past when adding e.g. a second kind of algebra solver

[avehtari]

If we would add integrate_1d_double_exponential I would then also add abs_tol argument and it would not be just an alias.

[WardBrian]

That would also be good, then.

If I recall correctly these functions are also already ready to be variadic at the math level but just aren’t in the language, should we make these new versions variadic from the start?

[avehtari]

I' m confident I can get Claude to create integrate_1d_double_exponential with added abs_tol argument, but I'm not confident I could trust it it to do the change to variadic arguments reliably. variadic arguments would be nice as there is a lot of packing and unpacking going on now, like this

    int K = x_i[1];
    real nu     = theta[1];
    real sigma  = theta[2];
    real sd1_1  = theta[3];
    real sd1_2  = theta[4];
    real rho    = theta[5];
    real z1     = theta[6 + 3 * K];

[WardBrian]

In the existing code (and this PR, looking over it briefly) the _impl function already is variadic. So we would just rename that to drop the _impl, remove the adaptor struct wrapping F, and we should be good to go.

I can help with this (and the stanc stuff, since it gets slightly more complicated to add a variadic), but I think it’s worth doing if we’re introducing new names anyway

[avehtari]

I can add integrate_1d_double_exponential with the additional arguments, and the test also if adding variadic arguments is easy. Should include them both to this PR?

[avehtari]

Ah, didn't see your latest comment. I can add integrate_1d_double_exponential with additional arguments and then you can add variadic

[WardBrian]

Sounds good. As long as claude follows the lead on the existing code, the extra steps to make it variadic should be very easy. The main reason it hasn’t been so far is that we need to come up with a new name for that version, which you’re already doing!