Observation time in ODE as parameters


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.

Comments and suggestions?


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


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).


Thanks, @betanalpha.