Help using new ode interface, getting incompatable types

Modeling
, ,
1

Hi,

Hopefully this is obvous to someone more familiar with the new ode interface.

I’ve worked through the examples provided by Upgrading to the new ODE interface but i’m still struggling to quite figure out what i’m doing wrong,
I’ve clearly got the types incorrect for the ode call, but I can’t see how?

I’ve simplifed a state-space model right down to what you see below.

I’m using cmdstan v2.26.1 / cmdstanr

the error is:

   -------------------------------------------------
    36:    for(i in 2:N){
    37:      // expects vector[] = ode_rk45(Func, vector y0, real t0, real[] times, )
    38:      C[i] = ode_rk45(predict, C[i-1], t0, dt[1], kw_obs, Csat);
                    ^
    39:    }
    40:  }
   -------------------------------------------------

Ill-typed arguments supplied to function 'ode_rk45'. Expected arguments:
(real, vector) => vector, vector, real, real[]
Instead supplied arguments of incompatible type:
(real, vector, real, real) => vector, real, real, real, vector, real

What the ode should be doing is predicting the value for C at each time-step, the time-step can be variable hence the length of dt.

functions {
  vector predict(real t,
                   vector y,
                   real kw,
                   real Csat){
    vector[1] dydt;
    dydt[1] = kw * (Csat - y[1]);
    return dydt;
  }
}

data {
  int<lower=2> N;
  vector[N] C_obs; // SML oxygen concentration (mmol m-3)
  vector[N] C_obs_error; // SML oxygen concentration error(mmol m-3)
  vector[N] kw_obs; // m s-1
  vector[N] Csat_obs; // SML Csat (mmol m-3)
  real dt[N];
}

transformed data{
  real t0 = 0;
}

parameters {
  real init_C;
  real Csat;
}

transformed parameters {
  vector[N] C; // latent state
  C[1] = init_C;
    // state model (no process error)
  for(i in 2:N){
    // expects vector[] = ode_rk45(Func, vector y0, real t0, real[] times, )
    C[i] = ode_rk45(predict, C[i-1], t0, dt[i-1], kw_obs, Csat);
  }
}

model {
    // priors
  Csat ~ normal(1, 1);
  init_C ~ normal(C_obs[1], C_obs_error[1]);
    // measurements
  C_obs ~ normal(C, C_obs_error);
}
2

I see three problems:

  • the state must be a vector, not a single number
  • solution times must be an array, not just a single number.
  • kw should be a single number but you’s passing in a vector.

Maybe this is what you want?

C[i] = ode_rk45(predict, [C[i-1]]', t0, {dt[i-1]}, kw_obs[i-1], Csat)[1,1];
3

On the other hand your ODE can be solved analytically as it has the canonical form du/dt=a + b * u.

4

Absolutely! I wish my full model was linear, that would be much easier.

2 Likes