Lotka-Volterra for three species

Modeling
, ,
1

Hi all, Stan novice here.

I’m trying to extend Bob Carpenter’s Lotka-Volterra case study to three species.

I’m however unable to make the sampling work and keep getting the following error:

SAMPLING FOR MODEL 'gLV3' NOW (CHAIN 1).
Unrecoverable error evaluating the log probability at the initial value.
Max number of iterations exceeded (1000).
[1] "Error in sampler$call_sampler(args_list[[i]]) : "
[2] "  Max number of iterations exceeded (1000)."

I’ve tried different initial values (below they are set to 0), different priors and different data sets but the error is always the same.

Any and all help is appreciated!

Stan code:


functions {
// function for the state equations
  real[] dz_dt(real t, 
               real[] z,  //state  
               real[] theta, //parameters
               real[] x_r, //data (unused) 
               int[] x_i) {
   real u = z[1];
   real v = z[2];
   real w = z[3];
   real r1 = theta[1];
   real r2 = theta[2];
   real r3 = theta[3];
   //
   real A11 = theta[4];
   real A12 = theta[5];
   real A13 = theta[6];
   real A21 = theta[7];
   real A22 = theta[8];
   real A23 = theta[9];
   real A31 = theta[10];
   real A32 = theta[11];
   real A33 = theta[12];
   
   real du_dt = (r1 + A11*u + A12*v + A13*w)*u;
   real dv_dt = (r2 + A21*u + A22*v + A23*w)*v;
   real dw_dt = (r3 + A31*u + A32*v + A33*w)*w;
   
   return { du_dt, dv_dt, dw_dt };
  }
}             

data {
  int<lower=0> N;      //number of measurements
  real ts[N];          //times
  real y0[3];          //initial measured populations
  real<lower=0> y[N, 3]; // measured populations
}

parameters {
  real theta[12];
  real<lower=0> z0[3];    //initial populations
  real<lower=0> sigma[3]; //measurement errors
}

transformed parameters {
  //population of measured years
  real z[N,3] 
    = integrate_ode_rk45(dz_dt, z0, 0, ts, theta,
                          rep_array(0.0, 0), rep_array(0,0), 
                          1e-6, 1e-5, 1e3);
}


model{
//priors
  sigma ~ lognormal(0, 0.5);
  theta[{1, 2, 3}] ~ normal(1, 1);
  theta[{4, 5, 6, 7, 8, 9, 10, 11, 12}] ~ normal(0, 0.5);

  z0[1] ~ lognormal(log(10), 5);
  z0[2] ~ lognormal(log(10), 5);
  z0[3] ~ lognormal(log(10), 5);
//likelihood
y0 ~ lognormal(log(z0), sigma);

for(k in 1:3) {
  y[ ,k] ~ lognormal(log(z[,k]), sigma[k]);
}
}

generated quantities {
  real y0_sim[3];
  real y_sim[N, 3];

  for (k in 1:3) {
    y0_sim[k] = lognormal_rng(log(z0[k]), sigma[k]);
    for (n in 1:N) {
      y_sim[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
    } 
  }
}

R code:


library(deSolve)
library(ggplot2)
library(rstan)

### Simulate data for gLV3
LotVmod <- function (Time, State, Pars) {
  with(as.list(c(State, Pars)), {

    ds1 = (r1 + A11*s1 + A12*s2 + A13*s3)*s1;
    ds2 = (r2 + A21*s1 + A22*s2 + A23*s3)*s2;
    ds3 = (r3 + A31*s1 + A32*s2 + A33*s3)*s3;
    
    return(list(c(ds1, ds2, ds3)))
  })
}

#Parameters:
Pars <- c(r1=0.6, r2=-0.8, r3=0.2,
          A11=0, A12=-0.03, A13=0,
          A21=0.03, A22=0, A23=0,
          A31=0, A32=0.3, A33=-0.2)
Time <- seq(0, 20, by = 1)
#initial state
State <- c(s1=10, s2=10, s3=10)
#make data frame
out <- as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))

#plot
# ggplot(data=out) + geom_line(aes(x=time, y=s1, colour="red")) + geom_line(aes(x=time, y=s2, colour="green")) + geom_line(aes(x=time, y=s3), colour="blue")


### STAN


# observation times
ts <- 1:(nrow(out)-1)
# number of observations
N <- nrow(out)-1 
# first observation
y0 <- c(out$s1[1], out$s2[1], out$s3[1])
# remaining observations
y <- as.matrix(out[2:nrow(out), 2:4])
# data
gLV3data <- list(N, ts, y0, y)

model <- stan_model("gLV3.stan")

inits <- function() list(r1=0, r2=0, r3=0,
                           A11=0, A12=0, A13=0,
                           A21=0, A22=0, A23=0,
                           A31=0, A32=0, A33=-0)

fit <- sampling(model, data=gLV3data, init=inits, seed=123)
2

Had a glance at this. Looks like something is going awry in the ODE integrate (it’s saying integrating between two of your time points took more than 1000 timesteps). If you throw in prints:

  print("print");
  z = integrate_ode_rk45(dz_dt, z0, 0, ts, theta,
                       rep_array(0.0, 0), rep_array(0,0),
                       1e-10, 1e-10, 1e4);
  print("sandwich");

You’ll see the ‘print’ but not the ‘sandwich’ on the output. Probly the easiest thing to do is print out the thetas and z0 and then try to integrate them in R and see what happens.

3

It looks like the max number of iterations was set to 1e4 (10,000), not 1e3 (1000), so I’m not sure what’s going on with that message. Are you using the default initializations? Did you simulate data from the model to use for testing?

I don’t know what would be causing that, I’m afraid.

Some unrelated comments that aren’t related to the problem you’re having:

vector[3] uvw = to_vector(z);
vector[3] r = to_vector(theta[1:3]);
matrix[3, 3] A = to_matrix(theta[4:12], 3, 3)';
vector[3] duvw_dt = (r + A * uvw) .* uvw;
return to_array_1d(duvw_dt);

and also

theta[1:3] ~ normal(1, 1);
theta[4:12] ~ normal(0, 0.5);
z0 ~ lognormal(log(10), 5);

Sorry there’s so much finagling of types required.

4

My bad, that was me that set it to 10k (just to test if a higher limit would work). @velait 's code has it at 1k (changing it is how I noticed the error)

5

Thank you for the comments. Sorry for the delayed response.

I’m not sure what happened but as of a few days ago the error doesn’t stop executing the program anymore.
There still are a few warnings in the beginning (see below) but sampling is done anyways.

Error evaluating the log probability at the initial value. Exception: lognormal_lpdf: Location parameter[1] is nan, but must be finite! (in 'model5ef375593d7f_gLV3' at line 70)

Sampling however takes considerable amounts of time on my laptop (1 chain and 300 iteration ~30min) and parameter estimates are far from correct values, which I guess is not surprising considering the low amount of iterations I used.

Even minor tweaking of one of the parameters (e.g A11 from 0 to -0.01 to 0.01) has drastic effects on the system so I wonder if the posterior distribution can have such a complex shape that in general sampling from it is not practically feasible.

6

Despite having written the case study, I’m not an expert on these models. I’m told they can be very chaotic in higher dimensions, so I’m not surprised they’re hard to fit.

Stan’s very sensitive to model fit. If models are badly misspecified for data, they can take a very long time to fit. Have you tried fitting data simulated directly from the model?

Can you send the program you’re using that’s slow? In what you have above, things like this are suboptimal:

y[ ,k] ~ lognormal(log(z[,k]), sigma[k]);

and can be helped a lot with some coaxing into better data structures. Just declaring y and z as matrices would be much faster than using 2D arrays, because matrices are column major indexed whereas arrays are implicitly row major. But presumably the ODE solver is consuming all the time, so that probably won’t help much.

One of the things that may be problematic here compared to the simple model where predator and prey are clearly defined is that the signs are uncertain in everything. Given that population sizes are positive, this can lead to uncertainty about which population is prey and which is predator w.r.t. the other species.

7

Hmm I see, I will have to read more about the topic.
The program I’m using is essentially the one in my first message, below with the suggested improvements. The data is simulated directly from the mode. Specifying initializations close to the correct values has no effect.

functions {
  // function for the state equations
  real[] dz_dt(real t, 
               real[] z,  //state  
               real[] theta, //parameters
               real[] x_r, //data (unused) 
               int[] x_i) {
                 
    vector[3] uvw = to_vector(z);
    vector[3] r = to_vector(theta[1:3]);
    matrix[3, 3] A = to_matrix(theta[4:12], 3, 3);
      
    vector[3] duvw_dt = (r + A * uvw) .* uvw;
      
    return to_array_1d(duvw_dt);
  }
}             

data {
  int<lower=0> N;      //number of measurements
  real ts[N];          //times
  real y0[3];          //initial measured populations
  matrix<lower=0>[N,3] y; // measured populations
  
}

parameters {
  real theta[12];
  real<lower=0> z0[3];    //initial populations
  real<lower=0> sigma[3]; //measurement errors
}

transformed parameters {
  //population of measured years
  real<lower=0>[N,3] z 
  = integrate_ode_rk45(dz_dt, z0, 0, ts, theta,
                       rep_array(0.0, 0), rep_array(0,0), 
                       1e-6, 1e-5, 1e4);
  //print(theta);
}


model{
  //priors
  sigma ~ lognormal(0, 0.5);
  theta[1:3] ~ normal(1, 1);
  theta[4:12] ~ normal(0, 2);
  
  z0 ~ lognormal(log(10), 5);
  
  //likelihood
  y0 ~ lognormal(log(z0), sigma);
  
  for(k in 1:3) {
    y[ ,k] ~ lognormal(log(z[,k]), sigma[k]);
  }
}

generated quantities {
  real y0_sim[3];
  real y_sim[N, 3];
  
  for (k in 1:3) {
    y0_sim[k] = lognormal_rng(log(z0[k]), sigma[k]);
    for (n in 1:N) {
      y_sim[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
    } 
  }
}

[edit: escaped code]

`