Parametrised integrals

Is it possible to define a function in stan that returns a parametrized integral? For example, an integrand parametrized by t_i,

I_i = \int_{a}^{b} f(s,t_{i}) \,\mathrm{d}s \,.

As far as I understand this is not possible, since the integrand has a strict definition

real integrand(real x, real xc, real[] theta,real[] x_r, int[] x_i)

that does not accept parametrized functions. Or the parametrization of the integrand should be fixed and dependent only on x_r and x_i.

My questions are:

  • Is my understanding correct?
  • Is there a workaround? I guess one workaround is to solve the integral using an ode solver.
  • Is there a particular reason the integrand has a strict definition? Since in the definition of the ode systems the user is free to pass any parameters.

Thank you for your time and your assistance.

That is what the real[] theta argument to the integrand is for.

Isn’t this variable only meant for variables of the block parameters?

It is not clear from my post, but in my case, t_i is not a parameter to be inferred, it is one of the data block.

If t is a real number, then pass it to the integrand function via x_r. If t is an integer, then use x_i.

I don’t think this can work since the integral is inside a loop and depends on t_i. I can pass the whole vector t through x_r but not the specific value t_i that I need on the i-th iteration.

Can’t you pass i through x_i, pass the whole vector \mathbf{t} through x_r, and pick out real t_i = x_r[x_i[1]]?

I tried this. I tried to change the value of x_i[1] inside the loop, i.e., x_i[1]=i and I got an error.

Is it allowed to change the value of transformed data inside the model block? If so, maybe I misinterpreted the error.

Integer arrays should always behave like data (the relevant distinction is not the block it’s defined in but whether the sampler needs to calculate a derivative). Something like this

real integrand(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
  real t_i = x_r[x_i[1]];
  return f(x, ti);
}
...
real x_r[N];
...
real theta[0];
for (i in 1:N) {
  ... = integrate_1d(integrand, a, b, theta, x_r, { i });
}

Note the curly braces that create a new array. Actually, this also works

integrate_1d(integrand, a, b, theta, {x_r[i]}, {1})

because anything built out of data-only expressions is still “data-like”.

1 Like

Hi Niko. Thanks for the reply. I tried it and I get the error:

   -------------------------------------------------
    51:        s[i] += 10.^gn * (x2[i]-1/x[i]);
    52:
    53:        real q = integrate_1d( J1_I, 0, pi(), {b,lw,th}, {x[i]}, {1} );
                        ^
    54:
    55:        s[i] += cm * q;
   -------------------------------------------------

Ill-typed arguments supplied to function 'integrate_1d'. Available signatures:
((real, real, real[], data real[], data int[]) => real, real, real, real[], data real[], data int[]) => real
((real, real, real[], data real[], data int[]) => real, real, real, real[], data real[], data int[], data real) => real
Instead supplied arguments of incompatible type: (real, real, real[], real[], int[]) => real, int, real, real[], real[], int[].

Same error with the integrate_1d(integrand, a, b, theta, x_r, { i }); approach.

From the error message, I understand that x_r a x_i must be data.

That’s weird because this compiles for me.

functions {
  real integrand(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
    real t_i = x_r[x_i[1]];
    return x + t_i;
  }
}
data {
  int N;
  real x_r[N];
  int x_i[0];
}
parameters {
  real b;
  real lw;
  real th;
}
model {
  for (i in 1 : N) {
    target += integrate_1d(integrand, 0, pi(), {b, lw, th}, {x_r[i]}, {1});
  }
}

Is the integrate_1d call inside a function? If so you must mark x as data in the signature.

Anyway, theta argument does not have to be from the parameters block.

integrate_1d( J1_I, 0, pi(), {b,lw,th,x[i]}, {1.0}, {1} )

Thank you Niko. This worked.

Yes, the integrate_1d is called inside another function and I have to add data in front of all the variables I entered in place of x_r in the integrate_1d call.