Bayesian inference to update a user-defined function from the data

Modeling
1

Inspired by Bob Carpenter’s blog on the Predator-prey population dynamics, I am trying to solve an inverse problem by inferring parameters from the data. A user defined parametric equation is denoted as y_r = a - b * exp(-y_e / c) + d * exp(-y_e / e). This is specified in the Stan code. a, b, c, d, and e are the parameters. A data set is labeled as y_r and y_e. But, I am stuck in errors. Can I get some help to address this problem?

// Stan code for the constitutive.R
functions {
  real[] cons_fn(real[] t,     // time
		 real[] theta,     // parameters in cons_fn
		 real y_r,
		 real y_e
		) 
  {
    // parameters
    real a = theta[1];
    real b = theta[2];
    real c = theta[3];
    real d = theta[4];
    real e = theta[5];

    // constitutive equation 
    real y_r  = a - b * exp(-y_e / c) + d * exp(-y_e / e);
    return {y_r, y_e};
  }
}

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

parameters {
  real<lower = 0> theta[5];   // {a, b, c, d, e }
  real<lower = 0> sigma[2];   // error scale
}

model {
  theta[{1, 2, 4}] ~ normal(1, 0.5, 0.5); // a, b, d
  theta[{3, 5}] ~ normal(0.05, 0.05); // c, e
  sigma ~ lognormal(-1, 1);
  y_init ~ lognormal(log(10), 1);
  for (k in 1:2) {
    y_init[k] ~ lognormal(log(y_init[k]), sigma[k]);
    y[k] ~ lognormal(log(z[k]), sigma[k]);
    //y[ , k] ~ lognormal(log(z[, k]), sigma[k]);
  }
}

generated quantities {
  real y_init_rep[2];
  real y_rep[N, 2];
  for (k in 1:2) {
    y_init_rep[k] = lognormal_rng(log(y_init[k]), sigma[k]);
    for (n in 1:N)
      y_rep[n, k] = lognormal_rng(log(z[n, k]), sigma[k]);
  }
}