Observation time in ODE as parameters

Developers
1

This is the first draft of allowing ts in ODE solvers to be parameters. The motivation is to make inference on the time of observations. What I did is to attach the RHS of the ODE to the decouple_states methods as that’s the sensitivity of the ODE variables with respect to the time.

The test model below shows what the implementation tires to achieve

functions {
  real[] sho(real t,
             real[] y, 
             real[] theta,
             real[] x,
             int[] x_int) {
    real dydt[2];
    dydt[1] = y[2];
    dydt[2] = -y[1] - theta[1] * y[2];
    return dydt;
  }
}
data {
  int<lower=1> T;
  real y0_d[2];
  real t0;
  real theta_d[1];
  real y_hat[T];
  real x[0];
  int x_int[0];
}
parameters {
  real<lower=0, upper=0.2> ts1;
  real<lower=ts1, upper=0.3> ts2;
  real<lower=ts2, upper=0.5> ts3;
}
transformed parameters{
  real<lower = 0> ts[T];
  for(i in 1:T) ts[i] = 0.1*i;
  ts[1] = ts1;
  ts[2] = ts2;
  ts[3] = ts3;
}

model {
  real sol[T,2];
  sol = integrate_ode_rk45(sho, y0_d, t0, ts, theta_d, x, x_int);
  y_hat ~ normal(sol[1:T, 2], 0.01);
}

With the ODE solution data, parameters ts1, ts2, and ts3 are inference results of the first three time steps.

ts1
ts1.png2099×2099 111 KB

ts2
ts2.jpg2099×2099 82.2 KB

ts3
ts3.jpg2099×2099 88.1 KB

Comments and suggestions?

2

It’d be good to get feedback from @charlesm93, @betanalpha, @billg, and @wds15. I don’t know enough about ODEs to be useful here.

3

Looks reasonable to me after going through the files, although it’s been a while since I’ve gone through them so take that with a grain of salt. My only comment is that in math/stan/math/prim/arr/functor/coupled_ode_system.hpp, Line 152 I recommend writing f(u(ti), ti) instead of f(ui, ti).

4

Thanks, @betanalpha.