Lotka-Volterra in LibBi

LV.bi

model LV {

  const h         = 0.02;   // time step
  const delta_abs = 1.0e-2; // absolute error tolerance
  const delta_rel = 1.0e-5; // relative error tolerance
  const epsilon   = 2.0e-2; // observation error tolerance

  dim s(2); // number of species
  dim n(2); // number of parameters for each species process

  param mu[s,n];    // mean of birth and death rates
  param sigma[s,n]; // volatility of birth and death rates
  noise w[s,n];     // noise terms
  state beta[s,n];  // birth and death rates
  state H, L;       // hares, lynxes
  obs y_H, y_L;     // observations

  sub parameter {
    mu[s,n] ~ uniform(0.0, 1.0);
    sigma[s,n] ~ uniform(0.0, 0.5);
  }

  sub proposal_parameter {
    mu[s,n] ~ truncated_gaussian(mu[s,n], 0.02, 0.0, 1.0);
    sigma[s,n] ~ truncated_gaussian(sigma[s,n], 0.01, 0.0, 0.5);
  }

  // The data are generated with parameters 0.5, 0.02, 0.4, 0.01 and
  // starting populations of 50 hares and 50 lynxes.
  sub initial {
    H ~ log_normal(log(50.0), 0.5);
    L ~ log_normal(log(50.0), 0.5);
    beta[s,n] ~ gaussian(mu[s,n], sigma[s,n]);
  }

  sub transition(delta = h) {
    w[s,n] ~ normal(0.0, sqrt(h));
    ode(h = h, atoler = delta_abs, rtoler = delta_rel, alg = 'RK4(3)') {
      dH/dt = beta[1,1] * H - beta[1,2] * H * L;
      dL/dt = beta[2,2] * H * L - beta[2,1] * L;
      dbeta[s,n]/dt = sigma[s,n] * w[s,n] / h;
    }
  }

  sub observation {
    y_H ~ log_normal(log(H), 0.2);
    y_L ~ log_normal(log(L), 0.2);
  }
}