LV.R

HaresLynxObservations <- structure(
  list(Year = 1900:1920, 
       Hares.x.1000 = c(30, 47.2, 70.2, 77.4, 
                        36.3, 20.6, 18.1, 21.4, 
                        22, 25.4, 27.1, 40.3, 57, 
                        76.6, 52.3, 19.5, 11.2, 
                        7.6, 14.6, 16.2, 24.7), 
       Lynx.x.1000 = c(4, 6.1, 9.8, 35.2, 59.4, 
                       41.7, 19, 13, 8.3, 9.1, 7.4, 
                       8, 12.3, 19.5, 45.7, 51.1, 
                       29.7, 15.8, 9.7, 10.1, 8.6)), 
  .Names = c("Year", "Hares.x.1000", "Lynx.x.1000"), 
  class = "data.frame", row.names = c(NA, -21L))

library(lattice)
xyplot(Hares.x.1000 + Lynx.x.1000 ~ Year, 
       data=HaresLynxObservations,
       auto.key=TRUE, t='l')

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores()) 
odeDso <- stan_model(file = "myLV_ODE.stan")

n <- nrow(HaresLynxObservations)

odeFit <- sampling(odeDso, 
                   data=list(T=n,
                             t0=-1/n,
                             y=HaresLynxObservations[, 2:3], 
                             ts=seq(0, 100, length=n)),
                   chains=1)

samples <- extract(odeFit, 'y_hat')

Y <- apply(samples[['y_hat']], c(2,3), mean)
plot(Y[,2] ~ Y[, 1], t="l")

myLV_ODE.stan

functions {
  real[] lv(real t,
    real[] y,
    real[] param,
    real[] x_r,
    int[] x_i) {
        real dydt[2];
        dydt[1] <-  param[1] * y[1] - param[2] * y[1] * y[2];
        dydt[2] <- -param[3] * y[2] + param[4] * y[1] * y[2];
        return dydt;
      }
}
data {
  int<lower=1> T;
  real y[T,2];
  real t0;
  real ts[T];
}
transformed data {
  real x_r[0];
  int x_i[0];
}
parameters {
  real y0[2];
  vector<lower=0>[2] sigma;
  real<lower=0, upper=1> param[4];
}
model {
  real y_hat[T,2];
  sigma ~ cauchy(0, 2.5);
  
  y0[1] ~ normal(0.5, 1);
  y0[2] ~ normal(1, 1);
  
  y_hat <- integrate_ode(lv, y0, t0, ts, param, x_r, x_i);
  for (t in 1:T){
    if (param[2] > 0.1)
    reject("param[2] too large: ", param[2]);
    if (param[4] > 0.1)
      reject("param[4] too large: ", param[4]);
    y[t] ~ normal(y_hat[t], sigma);
  }
}