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[1]; // Do not change
df_dx[1] = exp(-y^2 - x[1]^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[0];
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[2];
y_bounds[1] = ymin;
y_bounds[2] = 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?