Get ODE solution at specific time(s)

Modeling
1

Hi, I am working with an ODE. model I have data which is measured at specific time points of the time span. But I will need to solve the ODE for the entire time span and than need extract ODE output at observed time points for modelling. I would like to know how to change my code so that I could extract solutions of ODE at specific time span.
As an example, I use the following ODE model from Stan’s manual.

functions {
  real[] sho(real t,       //time
           real[] y,     // state
           real[] theta, //parameters
           real[] x_r,   // data (real)
           int[] x_i){  //data(integer)

     real dydt[2];

    dydt[1] = y[2];
    dydt[2] = -y[1] -theta[1] * y[2];

    return dydt;
   }
}
data {
  int<lower=1> T;
  real y[T, 2];
  real y0[2];
  real t0;
  real ts[T];
}
transformed data {
  real x_r[0];
  int x_int[0];
}
parameters {
  real theta[1];
  real<lower = 0> sigma;
}
transformed parameters{
real y_hat[T, 2];
  y_hat = integrate_ode_bdf(sho,y0, t0, ts, theta, x_r, x_int);
}
model {
   for (t in 1:T)
   y[t] ~ normal(y_hat[t], sigma);
   sigma ~ normal(0, 1);
   theta ~ normal(0, 1);
}

I use the following R code to simulate fake data.

sim_data = list(T = 10, y0 = c(1, 0), t0 = 0,
                ts = 1:10, theta = array(0.15,1))

fake_data = stan("sho.stan", data = sim_data,
            chains = 1, iter = 1,
            algorithm = "Fixed_param", seed =4838282)
y<- extract(fake_data, pars = "y_hat")$y_hat[1,,]

stan_data = list(T = 10, y0 = c(1, 0), t0 = 0,
                ts = 1:10, y =y )
fit <- stan("sho.stan", data =stan_data,
            chains = 4, iter = 1000,
            seed =4838282)

Now, suppose I have observed data at times 5:10, how would the above code be changed to model these observed data. The data block would be like:

  int<lower=1> T;
  real y[T, 2];
  real y0[2];
  real t0;
  real ts[T];
 real t_obs[5, 2]
real y_obs[5, 2]
}

and the data in the R code would be like:

stan_data = list(T = 10, y0 = c(1, 0), t0 = 0,
                ts = 1:10, y =y , y_obs = y[5:10, ], t_obs =data.frame(5:10, 5:10))

How should I modify the transformed data block so that I could extract ODE solution at the observed times . How should the model block be changed such that it models y_obs instead of y with mean as ODE solution corresponding to observed data.
THanks

2

Use T=5 and ts = 5:10, as ts is the time points where the ODE solution is retrieved into y_hat. ts is not the time step, which is determined adaptively.

3

Thanks @yizhang.