Stan code to fit Lotka-Volterra competition model

Lotka-Volterra_competition.stan

functions {
  vector dz_dt(real t, vector x, array[] real theta) {
    array[2] real r = theta[1:2];
    array[2] real K = theta[3:4];
    array[2] real alpha = theta[5:6];
    real dx1_dt = r[1] * x[1] * (1 - (x[1] + alpha[1] * x[2]) / K[1]);
    real dx2_dt = r[2] * x[2] * (1 - (x[2] + alpha[2] * x[1]) / K[2]);
    return [ dx1_dt, dx2_dt ]';
  }
}

data  {
  int<lower=0> N; // number of measurement times
  array[N] real ts; // measurement times > 0
  array[2] real<lower=0> y_init; // initial measured populations
  array[N, 2] real<lower=0> y; // measured populations
}

parameters {
  array[2] real<lower=0> r;
  array[2] real<lower=0> K;
  array[2] real<lower=0> alpha;
  vector<lower=0>[2] z_init; // initial population
  array[2] real<lower=0> sigma; // measurement errors
}

transformed parameters {
  array[6] real theta;
  theta[1:2] = r;
  theta[3:4] = K;
  theta[5:6] = alpha;
  array[N] vector<lower=0>[2] z = ode_rk45(dz_dt, z_init, 0, ts, theta);
}

model {
  r ~ normal(0, 2.5);
  K ~ normal(10, 5);
  alpha ~ normal(1, 1);
  sigma ~ normal(0, 2.5);
  z_init ~ lognormal(1, 2);
  for (k in 1:2) {
    y_init[k] ~ lognormal(log(z_init[k]), sigma[k]);
    y[ : , k] ~ lognormal(log(z[ : , k]), sigma[k]);
  }
}