1// Read the accompanying case study:
2// https://mc-stan.org/learn-stan/case-studies/lotka-volterra-predator-prey.html
3functions {
4 vector dz_dt(real t, // time
5 vector z, real alpha, real beta, real gamma, real delta) {
6 real u = z[1];
7 real v = z[2];
9 real du_dt = (alpha - beta * v) * u;
10 real dv_dt = (-gamma + delta * u) * v;
12 return [du_dt, dv_dt]';
13 }
14}
15data {
16 int<lower=0> N; // number of measurement times
17 array[N] real ts; // measurement times > 0
18 array[2] real y_init; // initial measured populations
19 array[N, 2] real<lower=0> y; // measured populations
20}
21transformed data {
22 real rel_tol = 1e-5;
23 real abs_tol = 1e-3;
24 int max_num_steps = to_int(5e2);
25}
26parameters {
27 real<lower=0> alpha, beta, gamma, delta;
28 vector<lower=0>[2] z_init; // initial population
29 array[2] real<lower=0> sigma; // measurement errors
30}
31transformed parameters {
32 array[N] vector[2] z = ode_rk45_tol(dz_dt, z_init, 0.0, ts, rel_tol,
33 abs_tol, max_num_steps, alpha, beta,
34 gamma, delta);
35}
36model {
37 alpha ~ normal(1, 0.5);
38 gamma ~ normal(1, 0.5);
39 beta ~ normal(0.05, 0.05);
40 delta ~ normal(0.05, 0.05);
41 sigma ~ lognormal(-1, 1);
42 z_init ~ lognormal(log(10), 1);
43 for (k in 1 : 2) {
44 y_init[k] ~ lognormal(log(z_init[k]), sigma[k]);
45 y[ : , k] ~ lognormal(log(z[ : , k]), sigma[k]);
46 }
47}
48generated quantities {
49 array[2] real y_init_rep;
50 array[N, 2] real y_rep;
51 for (k in 1 : 2) {
52 y_init_rep[k] = lognormal_rng(log(z_init[k]), sigma[k]);
53 for (n in 1 : N) {
54 y_rep[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
55 }
56 }
57}