Hi !

I’m currently working on a quite simple model for stan. It consists of time-series data and an ODE + noise model.

The equation for the ode relies on an input (the amount of “inducer”). In the current code (given below), the input is constant and given as data. What I would like to do is to give a time-varying input (also as data), for example two arrays, each one for the amount of one of the inducers at each time point, and solve the ODEs with it. Note that the amount of inducers is a step function After searching, i found several issues/post related to either using inference functions or the floor function to get inducers values from an array, but none actually gave a solution. How would you do this ?

**The code**

[details=Click to show the code] functions {

real hill_neg(real x,real theta,real eta) {

real aux=theta^eta;

return aux/(aux+x^eta);

}

real[] ode(real t,

real[] y,

real[] param,

real[] x,

int[] x_i) {

real dydt[2];

dydt[1] = hill_neg(y[2]*(1-x[1]),param[1],param[2])-y[1];
dydt[2] = hill_neg(y[1]*(1-x[2]),param[3],param[4])-y[2];

return dydt;

}

}

data {

int<lower=1> T;

real y[T,2];

real t0;

real ts[T];

real inducer[2];

}

transformed data {

int x_i[0];

}

parameters {

real<lower=0.0,upper=1> theta[2];

real<lower=1,upper=4> eta[2];

real<lower=0.00,upper=1> sigma;

}

transformed parameters {

real params[4] = {theta[1],eta[1],theta[2],eta [2]};

real vrnce = sigma^2;

}

model {

real y_hat[T,2];

//print("target (before) = ", target());

y_hat = integrate_ode_rk45(ode, {0.0,0.0}, t0, ts, params, inducer, x_i);

//print("y_hat = ",y_hat);

for (t in 1:T) {

y[t] ~ normal(y_hat[t], vrnce);

}

//print("target (after) = ",target());

}[/details]

Thanks a lot for your help