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

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]);
}
}
``````