Solving integrals that require solving an algebraic equation

Hello everyone, I have a purely mathematical problem which I wish to solve using Stan’s integrated tools but I am unable to do so.
Although this problem is not directly related to Stan, I am using it for the statistical modelling and I would like to keep it that way, hence the need of solving it using the tools it provides.

The notation used for this text will be the following:

  • \theta_i represents a parameter;
  • \phi_j represents a quantity that only depends on the parameters (i.e. \phi_j \equiv \phi_j (\theta) );
  • The functions f(t | \theta), g(t | \theta) and h(t | \theta) depend both on the time and the parameters.

In order to compute the quantities that can be directly compared with observations, I am required to find the value for the following integrals:
I_1 = \int_0^{\phi_2} \frac{1}{h(x | \theta)} dx
I_2 = \left( \int_{\phi_2}^\infty \frac{1}{h(x | \theta) f(x | \theta) } dx \right)

Where h(t, \theta) can be obtained by numerically solving the following equation:
(h^2 - 2 \phi_1)e^{\phi_1/h^2} = g(t | \theta)               (1)

Where both f(t | \theta) and g(t | \theta) are fully determined and analytical (not shown here for brevity).

I’ve asked a similar question before in “Unable to integrate function obtained from algebraic solver”, where a user suggested me to use the Stan’s DAE solver, which works very well to obtain the value of I_1.
All you have to do is re-write this system of equations as a system of differential algebraic equations and perform a variable change because Stan’s DAE solver doesn’t allow time as parameters, as clarified in “Compilation error on using the Differential-Algebraic equation (DAE) solver on time that depends on the parameters” thanks once again to the wonderful support provided by users of this forum.
Doing so results in the following system of equations:
(y_1^2 - 2 \phi_1)e^{\phi_1/y_1^2} - g(t | \theta) = 0
\frac{dy_2}{dt} - \frac{\phi_2}{y_1} = 0

Where \vec{y} is the quantity we are solving for and its components are defined by
y_1 (t | \theta) \equiv (h(t|\theta)
y_2 (t | \theta) \equiv \int_0^{\phi_2} \frac{1}{h(x|\theta)}dx

Using the initial conditions:
y_1 (0 | \theta) = 1
\frac{y_1}{dt} (0 | \theta) = \frac{1}{2e^{\phi_1}} \frac{3(\theta_2 + \theta_3) + 4\theta_4}{1 - \phi_1 + 2\phi_2^2}
y_2 (0 | \theta) = 0
\frac{y_1}{dt} (0 | \theta) = 1

We then call Stan’s built in DAE solver to solve this system from t = 0 to t = 1.

If anybody is interested in looking at the full model file (which is running as expected, so it’s not really a requirement I’m just sharing it here just because), but with the names of the parameters/variables different from those used here, can do so in the following file:
model.stan (2.7 KB)

However, I don’t think I’m able to re-write I_2 in such a way that I can use this approach, using the DAE solver, because one of its bounds is infinity instead of 0.

An approach I tried to take was to re-write equation (1) into an ODE and then use the ode_rk45 to obtain numerically h(t|\theta) and then obtain I_1 and I_2 using integrate_1d, but unfortunately I came across a segmentation fault and reported it in Segmentation fault on making use of `ode_rk45` with `integrate_1d` · Issue #2848 · stan-dev/math · GitHub, with no solution found so far

This leaves me with my question:

  • How can I numerically find the value for both I_1 and I_2 using the tools provided by Stan?

Thanks for the careful question and math and links to the full model. This one’s a bit esoteric for most of our devs, including me. I think @charlesm93 may be able to help if he sees this ping.

Many thanks for your reply.
Earlier today I found a viable solution.

The trick I used is to use an approximate form of equation (1) instead of the full blown numerical solution.

Because the value of \phi_2 is, for a reasonable region in parameter space, a very very high value, something which I confirmed independently using Mathematica, this means that the integrand of I_2 will be computed for very high values.
This corresponds to the case where h is very large, meaning that the factor \phi_1/h^2 in equation (1) is very small.
Performing a Taylor series expansion to first order of the previous factor I was able to get an analytical solution for h and using that to compute I_2, while using the full blown numerical solution to compute I_1 to avoid errors due to approximations.

Now, in my attempt of simplifying the problem I didn’t give a few details which in this case were necessary to suggest such a solution.
Maybe there is another way of solving this without relying on approximate solutions, but all my analysis seems to indicate that the error is completely negligible.

I can only thank you and the other members of the Stan development team even bothering to read random questions asked by random people, which in this case aren’t even related to Stan!
It’s a great tool and this forum is a great place as well!