Yeah, integrate_1d can be rather frustrating.
The integrate_1d function is basically a wrapper around the Boost 1d quadrature stuff over here: https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/quadrature.html .
On top of doing the integrals, it will also compute gradients of these integrals which can blow up. But I don’t think it does that for expose_stan_functions.
To understand the error, check out: https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/double_exponential/de_tol.html
I’m not sure what algorithm R’s integrate uses.
Looking at the function, Phi can be a rather fragile function to work with. I wonder if that is the thing blowing up. Stan and R use different implementations of functions like this. Numeric overflows/underflows mess up the integrator. I guess the place to start would be checking if either the pdf or cdf in your integral are failing to evaluate at any point.
Hopefully it’s one of those functions messing up so we can rewrite the integral and it will work (like what happened here: Using integration to fit the difference of gamma distributed random variable)