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?