Solve ODE system with a known vector of one compartment

Borrowing the ODE example from the stan documentation, I have an ODE system with one know time-dependent vector (vector Z[T] = [integers]). I know I need an interpolation function to use Z_t[i] in each step of ODE solver and use in transformed parameters as well. I have changed the code as below but cannot define the function and make it work:


functions {
real[] intropolation(){
}
real[] sir(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i) {

real 
real dydt[3];

z_t= intropolatin()
dydt[1] = -theta[1]* y[2] * y[1] ;
dydt[2] = theta[1]* y[2] * y[1] - theta[2]* y[2] + z_t;
dydt[3] =  theta[2]*  y[2];

return dydt;
}
}
data {
int<lower=1> T;
real y[T,3];
real t0;
real ts[T];
vector[T] Z;
}
transformed data {
int x_i[0];
real x_r[T];
for (i in 1:T){
    x_r[i] = Z[i];
}
transformed parameters{
  real y[T, 3];
  real cases[T - 1];
  {
    y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
  }
  for (i in 1:T-1){
   cases[i] =  y[i, 1] - y[i+1, 1] + Z[i];
  }
}

Thank you in advance for your help.
#modeling #Interpolation #ODE solver #Bio

Hi,
unfortunately the question is a bit unclear, could you be more specific in what you mean by “cannot define the function and make it work” - do you get a syntax error? Are you unsure on what the function should actually do? Can you run the model, but the results look weird?

Additionally, in most cases like that, you usually don’t have the “true” values of the vector, but rather have observations of the true value, so it might be more sensible to treat the true values of z as an unknown spline (or Gaussian process, but that’s a bit more challenging) - which you can readily evaluate at arbitrary time points - and then formulate an observation model that constraints the spline to fit the observed value. Would that make sense in your situation or do you believe that the observation noise is indeed negligible for your use case?

Best of luck with your model.

2 Likes