Separating the outputs of ODE solver

Modeling
, , ,
1 
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.

2

See 5.6 Array data types | Stan Reference Manual about accessing elements of arrays. The output of the ode solver is an array of vectors.

Without knowing the data it’s hard to know why the two params are not identified correctly. A similar situation I ran into before is that y_1^2 << y_2 so the first two terms become unidentifiable. Obviously that’s only a guess.

1 Like
3

Thanks for the reply. If it helps, I can simply post the data here. There is no such consistent difference, positive or negative, between y_1^2 and y_2. I will take a look into the array documentation you linked, thanks for that.

And here is the data in CSV format. Any help in debugging is highly appreciated. Thank you so much.

t,y1,y2
0,2,0
1,1.50809694334058,-0.78030462089834
2,0.323175864421418,-1.83310134300764
3,-1.86840649338954,-1.0225792347164
4,-1.74357133277781,0.624519634716855
5,-0.840045152717157,1.30511332360852
6,1.27461859803541,2.44657186967907
7,1.92219626388251,-0.434789677505298
8,1.21698287932375,-0.984612109163882
9,-0.406238046128214,-2.52323554208559
10,-2.01193364713389,0.0301202670980852

4

I wonder the rationale of the constraints and uniform priors. Are the constraints based on actual knowledge of the dynamical system the ode describes? Are the uniform priors confirmed by the prior predictive checks?

If not maybe you can loose them, something like

  real<lower=0> alpha;
  real<lower=0> eta;
  real<lower=0> lambda;
  real<lower=0> beta;
//...
  alpha ~ normal(3, 1);
  eta ~ normal(0, 2);
  lambda ~ normal(0, 2);
  beta ~ normal(4, 2);
//...

Here I assume positive params and changed the prior based based on the uniform distr interval.

1 Like
5

Hmm, seems like a very good suggestion. Because now, the initial conditions are being estimated correctly. Although, all the parameters are wrong. This may be because this is really a toy example and does not really represent any real physical phenomenon.

Do you think even then the parameters should be estimated (at least to some degree) correctly?

I cannot thank you enough for your valuable inputs and help on this matter. Thanks.

6

I’m guessing what happened is that your hard-coded constraints have excluded the “true” initial conditions from the parameter space, while the new model no longer suffers from that it is unidentifiable with respect to \alpha,\beta,\eta, \lambda, and the priors I gave is biased toward the “wrong” ones. This is why I was asking about the dynamics behind the ODE(what physics/scientific phenomenon the ODE describes), because only with some understanding of the system can one enrich the model with informative priors and bias toward the “correct” param values.

1 Like
7

Thanks for all your valuable inputs on this. I have some more information on the system. Seems like this toy data was constructed from the Van der Pol equation, modifying the parameter values a little bit. What do you think might be going wrong with my code and estimation of the parameters in the light of this new information?

Thanks.