# Data vs reference parameters for ODE integrators

**URL:** <https://discourse.mc-stan.org/t/data-vs-reference-parameters-for-ode-integrators/9671>\
**Category:** Modeling\
**Created:** [July 11, 2019, 1:30pm UTC](https://discourse.mc-stan.org/t/data-vs-reference-parameters-for-ode-integrators/9671 "2019-07-11T13:30:05Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![NougatRillettes](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nougatrillettes/32/364_2.png) [@NougatRillettes](https://discourse.mc-stan.org/u/NougatRillettes)\
**Post date:** [July 11, 2019, 1:30pm UTC](https://discourse.mc-stan.org/t/data-vs-reference-parameters-for-ode-integrators/9671/1 "2019-07-11T13:30:05Z")

</div>

Hi everyone!

This post is somehow related to an old [post](https://discourse.mc-stan.org/t/time-varying-input-for-ode-solving-aka-forced-ode/989) of mine.

I need to fit experimental data to an ODE model with _time-varying piecewise continuous inputs_. I just dived into it and tried to implement a simple loop that would call the ODE solver between each of the time points, but I am having lots of issues figuring out how to code it without having issues with the inferred non-“dataness” of my time-points.

From what I gathered from this [message](https://discourse.mc-stan.org/t/ode-integrator-calls-from-user-defined-functions/2043/6) this is because local variables can never be inferred as “dataish”, but I can’t see how to code this.

The full code is below but it breaks down as follow:

## The data

Data is in the form of 4 vectors:

- `inputdurations` containing the duration of each step,
- `inputvalues` containing the value of each step
- `fluotimes` containing each measurement’s time
- `fluovalues` containing each measurement (fluorescence) value

## The model

It is a deterministic ODE model with additive Gaussian noise. There are 3 state variables (of which one is measurable) and 10 parameters.

## The implementation

> **The code**
>
> ```stan
> functions {
> real hillpos(real x, real theta, real eta) {
> return x^eta/(x^eta + theta^eta);
> }
> 
> real[] dydt(real t, real[] y, real[] prms, real[] x_r, int[] x_i) {
> real kon = prms[1];
> real koff = prms[2];
> real kb = prms[3];
> real kmax = prms[4];
> real kd = prms[5];
> real eta = prms[6];
> real kdegr = prms[7];
> real ktrans = prms[8];
> real kdegp = prms[9];
> real tftot = prms[10];
> 
> real u = x_r[1];
> 
> real tfon = y[1];
> real mrna = y[2];
> real prot = y[3];
> 
> real res[3];
> res[1] = u*kon*(tftot-tfon)-koff*tfon;
> res[2] = kb + kmax*hillpos(tfon,kd,eta)-kdegr*mrna;
> res[3] = ktrans*mrna-kdegp*prot;
> return res;
> }
> 
> real mul_gaussian_lpdf(real y, real mu, real sigma) {
> return normal_lpdf(y | mu, sigma*mu);
> }
> }
> 
> data {
> int<lower=0> ninput;
> vector[ninput] inputdurations;
> vector[ninput] inputvalues;
> 
> int<lower=0> nfluo;
> vector[nfluo] fluotimes; // /!\ Assumed to be increasing order!
> vector[nfluo] fluovalues;
> }
> 
> transformed data {
> real tbeg = 0;
> }
> 
> parameters {
> real<lower=0> kon;
> real<lower=0> koff;
> real<lower=0> kb;
> real<lower=0> kmax;
> real<lower=0> kd;
> real<lower=0> eta;
> real<lower=0> kdegr;
> real<lower=0> ktrans;
> real<lower=0> kdegp;
> real<lower=0> tftot;
> 
> real<lower=0> sigma;
> }
> 
> model {
> real y0[3] = {0.0, 0.0, 0.0};
> int nexttidx = 1;
> real params[10] = {kon, koff, kb, kmax, kd, eta, kdegr, ktrans, kdegp, tftot};
> int x_i[0];
> 
> kon ~ mul_gaussian(0.0016399, 0.1);
> koff ~ mul_gaussian(0.34393,0.1);
> kb ~ mul_gaussian(0.02612,0.1);
> kmax ~ mul_gaussian(13.588,0.1);
> kd ~ mul_gaussian(956.75,0.1);
> eta ~ mul_gaussian(4.203,0.1);
> kdegr ~ mul_gaussian(0.042116,0.1);
> ktrans ~ mul_gaussian(1.4514,0.1);
> kdegp ~ mul_gaussian(0.007,0.1);
> tftot ~ mul_gaussian(2000,0.1);
> 
> for (i in 1:ninput) {
> real val = inputvalues[i];
> real dur = inputdurations[i];
> real tend = tbeg + dur;
> real oderes[1,3];
> real taux = tbeg;
> 
> while ((nexttidx <= nfluo) && (fluotimes[nexttidx] <= tend)) {
> real t = fluotimes[nexttidx];
> real y = fluovalues[nexttidx];
> nexttidx += 1;
> 
> oderes = integrate_ode_rk45(dydt, y0, taux, {t}, params, {val}, x_i);
> y0 = oderes[1];
> taux = t;
> y ~ normal(oderes[1,3], sigma);
> }
> if (taux < tend) {
> oderes = integrate_ode_rk45(dydt, y0, taux, {tend}, params, {val}, x_i);
> y0 = oderes[1]
> }
> tbeg = tend;
> }
> }
> 
> ```

> **Error message**
>
> ```
> SYNTAX ERROR, MESSAGE(S) FROM PARSER:
> Third argument to integrate_ode_rk45 (initial times) must be data only and not reference parameters.
> Fourth argument to integrate_ode_rk45 (solution times) must be data only and not reference parameters.
> Sixth argument to integrate_ode_rk45 (real data) must be data only and not reference parameters.
> 
> ```

---

<div class="post-metadata">

**Author:** ![yizhang](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/yizhang/32/15673_2.png) [@yizhang](https://discourse.mc-stan.org/u/yizhang)\
**Post date:** [July 11, 2019, 4:14pm UTC](https://discourse.mc-stan.org/t/data-vs-reference-parameters-for-ode-integrators/9671/2 "2019-07-11T16:14:13Z")

</div>

Yes, everything declared in `model` block becomes a parameter, which is why Stan is complaining things declared there are fed to data-only args. To do time-varying ODE I’d modify the ODE function definition to reflect that the right-hand-side behaves differently in different time and/or parameter.