```
functions {
vector test(real t,
vector y,
real alpha,
real eta,
real lambda,
real beta) {
vector[2] dydt;
dydt[1] = alpha * y[2];
dydt[2] = eta * y[2] - lambda * y[1]^2 * y[2] - beta * y[1];
return dydt;
}
}
data {
int<lower=1> T;
array[T] vector[2] y;
real t0;
array[T] real<lower=1> ts;
// vector[2] y0;
}
parameters {
vector[2] y0;
vector<lower=0>[2] sigma;
real<lower=2, upper=4> alpha;
real<lower=0, upper=3> eta;
real<lower=0, upper=3> lambda;
real<lower=3, upper=6> beta;
}
model {
vector[2] mu[T] = ode_rk45(test, y0, t0, ts, alpha, eta, lambda, beta);
sigma ~ cauchy(0, 1);
alpha ~ uniform(2, 4);
eta ~ uniform(0, 3);
lambda ~ uniform(0, 3);
beta ~ uniform(3, 6);
y0 ~ normal(0, 1);
for (t in 1:T) {
y[t] ~ normal(mu[t], sigma);
}
}
```

This is the code that I have written for fitting an ODE given by the equations:

\frac{dy_1}{dt} = \alpha y_2 \\
\frac{dy_2}{dt} = \eta y_2 - \lambda y_1^2 y_2 - \beta y_1

with the hyper-priors as the same as specified in the model block. Now, I am trying to fit the likelihoods as follows:

y_1^*(t) \sim N(y_1(t), \sigma_1 = (0.1*y_1(t))) \\
y_2^*(t) \sim N(y_2(t), \sigma_2 = (0.1*y_2(t)))

But, I could not find a way to separate the integrated values, viz., y_1(t) \ \& \ y_2(t) from my first line of the model block. As a result, I had to use a diffuse prior on the \sigma's, instead of the way I want them to be.

To give a bit more clarification and context to this, I know the true values of the parameters and the initial states, and I am trying to get estimates as close to the known truth as possible. That way, I will know that the model-fitting and estimation is working as expected with the Stan and R.

But this model run gives reasonable posterior estimates for the \alpha and the \beta but are way off from the reality in case of \eta, \lambda and the initial states. I tried providing the initial values as known, instead of estimating them, even then the 4 parameter values donâ€™t change much at all from the (given here) model estimates where I am estimating the initial states.

Is there any way that I can achieve this kind of fitting using Stan and R, where I can separate the \sigma's with the help of separating the derivatives and integrate them separately in the model block? Also, if there are any inconsistencies in my code and what I want to do, please point that out to me. Any help in the right direction will be hugely appreciated! Thank you very much.