# Passing data to the ODE integrator in user-defined functions for double integration

I have written user-defined functions to perform the double definite integration \int_{x_0}^{x_1}\int_{y_0}^{y_1}f(x, y)dydx. Although it works, the bounds of the integration have to be written into the function body in unnatural ways.

functions {
// Example integrand
real[] integrand_ode(real y, real[] f, real[] theta, real[] x, int[] x_i) {
real df_dx; // Do not change
df_dx = exp(-y^2 - x^2); // Function of x and y to be integrated
return(df_dx);
}
//User-defined function of x after y integrated out
real[] first_integral(real x, real[] f, real[] theta, real[] x_r, int[] x_i) {
return(integrate_ode_rk45(integrand_ode,
rep_array(0.0, 1),
0.0, // ymin
rep_array(1.0, 1),
rep_array(0.0, 0),
rep_array(x, 1),
x_i)[1,]);
}
// Final integration wrt x
real second_integral(real xmin, real xmax ) {
int x_i;
return(integrate_ode_rk45(  first_integral,
rep_array(0.0, 1),
xmin,
rep_array(xmax, 1),
rep_array(0.0, 0),
rep_array(0.0, 0) ,
x_i)[1,1]);

}

}
data {
}
model {}



What I want is for second_integral to look something like this

 real second_integral(real xmin, real xmax , real ymin, real ymax) {

real y_bounds;
y_bounds = ymin;
y_bounds = ymax;

return(integrate_ode_rk45(  first_integral,
rep_array(0.0, 1),
xmin,
rep_array(xmax, 1),
rep_array(0.0, 0),
y_bounds ,
rep_array(0, 0))[1,1]);

}


so that the bounds are all entered in one place and transmitted as data to the functions. But if I parse that form of the function I get:

sixth argument to integrate_ode_rk45 (real data) must be data only and not reference parameters error

That error message has been dealt with in this forum https://discourse.mc-stan.org/t/standalone-use-of-ode-integrator-in-r/1972 but only in the context that the sixth argument is empty.

Is it possible to do what I want?

Very q: Probably it is possible to do, yes… BUT: all temporaries inside user functions of stan are casted to being parameters. So even the rep_array(0,0)[1,1] is for Stan a parameter and not data.

The solution is to define all your data arguments inside the data or the transformed data and then pass these into the functions. You are allowed to subset these, but not more. Then it will stay as data and it should work.

Thanks for such a quick response. You will have noticed that my code draws heavily on your excellent single integrator function.
I am not quite clear what you are suggesting in this case. In my specific application for this function all the bounds of the integrals will be in the data, so that part is straightforward. What I don’t understand is what you mean by

pass these into the functions

I should be most grateful if you could be more specific.

before you get a double integration working, please get a single integration working which is called from inside a user function. I don’t have an immediate example which I can refer you to, but this more simpler exercise should let things sink in, I hope.

Two dimensional integration is going to get very expensive and relatively inaccurate, even for seemingly simple integrands. Effective two dimensional quadrature is more than just iterating effective one-dimensional quadratures, and there really aren’t effective three-dimensional quadrature methods at all!