Revamp Lotka-Volterra ODE model with new syntax

Modeling
1

Hi!
I am trying to translate the Lotka Volterra Predator-Prey population dynamics tutorial using the Stan data types.
But If I define

functions{
  vector lotka_volterra(real t, real[] z, real[] theta) {
    real u = z[1];
    real v = z[2];
    vector[2] dz_dt;

    real alpha = theta[1];
    real beta = theta[2];
    real gamma = theta[3];
    real delta = theta[4];

    dz_dt[1] = (alpha - beta * v) * u;
    dz_dt[2] = (-gamma + delta * u) * v;
    return dz_dt;
  }
}

and defining:

data {
  int<lower = 0> N;           // num measurements
  array[N] real ts;           // measurement times > 0
  vector[2] y_init;           // initial measured population
  array[N] vector[2] y;       // measured population at measurement times
}

parameters {
  vector[4] theta;   // theta = { alpha, beta, gamma, delta }
  vector[2] z_init;  // initial population
  vector<lower = 0>[2] sigma;   // error scale
}

and then use it in

transformed parameters {
  array[N] vector[2] z = integrate_ode_rk45(lotka_volterra, z_init, 0, ts, theta);
}

I get the following compilation error:

Ill-typed arguments supplied to function 'integrate_ode_rk45':
(<F1>, vector, int, array[] real, vector)
where F1 = (real, vector, vector) => vector
Available signatures:
(<F2>, array[] real, real, array[] real, array[] real, data array[] real,
 data array[] int, data real, data real, data real) => array[,] real
where F2 = (real, array[] real, array[] real, data array[] real,
            data array[] int) => array[] real
  Expected 10 arguments but found 5 arguments.
(<F2>, array[] real, real, array[] real, array[] real, data array[] real,
 data array[] int) => array[,] real
  Expected 7 arguments but found 5 arguments.

Where’s the error?
Thanks in advance…
BTW I tag @Bob_Carpenter as the original author of the case study

2 Likes
2

the ode rhs function seems to be missing the x_r and x_i arguments, I think. Have a look at the refs for this in the manual…and consider upgrading this to the variadic interface Upgrading to the new ODE interface … for that you need a more recent Stan though.

3

We have an issue now with most of the case studies not working with 2.33+:

1 Like
4

It would be great to see some general advice about this. As a package developer I’d like my code to be robust to the Stan versions that users have but 2.33 has caused some issues. I presume this is too tough for the autoformatter to handle?

5

I doubt it can upgrade ODE interfaces for you (but I honestly don’t know, worth checking). But it can handle syntax issues like the new array syntax, _log to _lpdf, # to //, etc.

1 Like
6

Actually x_r and x_i were the missing pieces and now the source code compiles. And I missed all the implications of the variadic interface. So, thanks @wsd15 for the link to the documentation!

I post here the updated model for future reference (@jonah: could the code be of any help in updating the case studies?) or for anyone else in need. I hope you’ll find it useful and a good starting point for you to work on it. Just pay attention to priors: they must be set accordingly to the problem you will solve!

functions{
  vector lotka_volterra(real t, vector z, real alpha, real beta, real gamma, real delta) {
    vector[2] dydt;

    dydt[1] = (alpha - beta * z[2]) * z[1];
    dydt[2] = (-gamma + delta * z[1]) * z[2];

    return dydt;
  }
}

data {
  int<lower = 0> N;           // num measurements
  array[N] real ts;           // measurement times > 0
  vector[2] y_init;           // initial measured population
  array[N] vector[2] y;       // measured population at measurement times
}

parameters {
  real alpha;
  real beta;
  real gamma;
  real delta;
  vector[2] z0;                 // initial population
  vector<lower = 0>[2] sigma;   // error scale
}

transformed parameters {
  array[N] vector[2] z = ode_bdf(lotka_volterra, z0, 0, ts, alpha, beta, gamma, delta);
}

model {
  alpha ~ std_normal();
  beta ~ std_normal();
  gamma ~ std_normal();
  delta ~ std_normal();
  sigma ~ std_normal();
  z0 ~ std_normal();

  for (k in 1:2) {
    y_init[k] ~ normal(z0[k], sigma[k]);
    y[ , k] ~ normal(z[, k], sigma[k]);
  }
}
2 Likes