# ODE with time-varying parameters whose trajectory is just polynomial

**URL:** <https://discourse.mc-stan.org/t/ode-with-time-varying-parameters-whose-trajectory-is-just-polynomial/33090>\
**Category:** Modeling\
**Tags:** rstan, specification\
**Created:** [October 18, 2023, 3:38pm UTC](https://discourse.mc-stan.org/t/ode-with-time-varying-parameters-whose-trajectory-is-just-polynomial/33090 "2023-10-18T15:38:57Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mbrabec](https://avatars.discourse-cdn.com/v4/letter/m/d6d6ee/32.png) [@mbrabec](https://discourse.mc-stan.org/u/mbrabec)\
**Post date:** [October 18, 2023, 3:38pm UTC](https://discourse.mc-stan.org/t/ode-with-time-varying-parameters-whose-trajectory-is-just-polynomial/33090/1 "2023-10-18T15:38:57Z")

</div>

I am trying to fit a system of two ODEs model with two time-varying parameters. Trajectory of these parameters is approximated by polynomials. Using rstan.

When compiling, I get the message:  
Error in stanc(file = file, model\_code = model\_code, model\_name = model\_name, :  
0

Syntax error in ‘string’, line 56, column 53 to column 54, parsing error:

Ill-formed expression. Found identifier. There are many ways to complete this to a well-formed expression.

When I tried to modify the code, simplify it to only one ODE with one time-varying parameter, the problem remains the same (only the line # shifts), so that it seems that I miss something deep and/or obvious. Could you point me somehow what it can be?  
I am listing my Stan code below.  
By the way, the line 56 is within the transformed parameters block and looks comletely innocent:  
vector \<lower=0\> [n\_radon] mu\_rn;  
This increases my confusion substantially (seems that the problem propagates from some other places, but could not figure from where) …

All the Best  
Marek B.

```stan
functions {
  real k(real t, vector beta_k){
    real k0 = beta_k[1];
    real k1 = beta_k[2];
    real k2 = beta_k [3];
    real k3 = beta_k[4];
    real k4 = beta_k [5];
    real k5 = beta_k [6];
    real k6 = beta_k[7];
    real k7 = beta_k [8];
    real k8 = beta_k [9];
    k=1/(1+exp(-(k0+k1*t+k2*t^2+k3*t^3+k4*t^4+k5*t^5+k6*t^6+k7*t^7+k8*t^8)));
    return k;
  }
  real G(real t, vector beta_G){
    real G0 = beta_G[1];
    real G1 = beta_G[2];
    real G2 = beta_G [3];
    real G3 = beta_G[4];
    real G4 = beta_G [5];
    real G5 = beta_G [6];
    real G6 = beta_G[7];
    real G7 = beta_G [8];
    real G8 = beta_G [9];
    G=exp(G0+G1*t+G2*t^2+G3*t^3+G4*t^4+G5*t^5+G6*t^6+G7*t^7+G8*t^8);
    return G;
  }
  vector sysode(real t, real y, vector beta_k, vector beta_G, real G_tracer) {  
    real tracer =y[1];
    real radon = y[2];
    vector[2] dy_dt;
    dy_dt[1] = G_tracer-k*tracer;
    dy_dt[2] = G-k*radon;
    return dy_dt;
  }
}
data {
  int < lower =0> n_tracer;
  int < lower =0> n_radon;
  vector [n_tracer] tracer;
  vector [n_radon] rn;
  vector [n_tracer] times_tracer;
}
parameters {
  real<lower=0> sigma_tracer;
  real<lower=0> sigma_rn;
  real<lower=0> G_tracer;
  vector [9] beta_k;
  vector [9] beta_G;
  vector [2] y0;
}
transformed parameters {
  vector <lower=0> [n_tracer] l_tracer;
  vector <lower=0> [n_tracer] l_rn;
  vector <lower=0> [n_radon] mu_rn;
  array[n_tracer] vector[2] y = ode_rk45(sysode, y0, t0=0, ts=times_tracer, beta_k, beta_G, G_tracer);
  l_tracer=y[,1];
  l_rn=y[,2];
  for(j in 2:(n_radon-1)){
    mu_rn[j]=sum(l_rn[((j-2)*12+2):((j-1)*12+1)]);
    //mu_rn[j]=sum(l_rn[((j-2)*12+2):((j-1)*12+1)])/12;
  }
  mu_rn[1]=mu_rn[2];
  mu_rn[n_radon]=sum(l_rn[((n_radon-2)*12+2):((n_radon-1)*12)]);
  //mu_rn[n_radon]=sum(l_rn[((n_radon-2)*12+2):((n_radon-1)*12)])/12;
}
model {
  G_tracer ~ normal(0, 10);
  beta_k~normal(0,1);
  beta_G~normal(0,10);
  sigma_tracer ~ normal(0, 10);
  sigma_rn ~ normal(0, 10);
  y0[1]~lognormal(log(50),log(50));
  y0[2]~lognormal(log(500),log(500));

  for(i in 1:n_tracer){
    tracer[i] ~ normal (l_tracer[i], sigma_tracer);
  }
  for(j in 1:n_radon){
    rn[j] ~ normal (mu_rn[j], sigma_rn);
  }
}

```

[edit: escaped Stan code and auto-indented]

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [November 1, 2023, 9:20pm UTC](https://discourse.mc-stan.org/t/ode-with-time-varying-parameters-whose-trajectory-is-just-polynomial/33090/2 "2023-11-01T21:20:53Z")

</div>

That is not a good error message from our compiler, but it is the right line:

```stan
  array[n_tracer] vector[2] y = ode_rk45(sysode, y0, t0=0, ts=times_tracer, beta_k, beta_G, G_tracer);

```

Stan doesn’t allow named arguments, so you can’t have `t0 = 0` or `ts = times_tracer`. Just remove the name and equal sign and this problem should go away.

I would also strongly recommend adding spaces to the code and not renaming all the elements of an array like in the definition of `real G`.

I’d also strongly recommend using built-ins for things like inverse logit, so rather than

```stan
k = 1 / (1 + exp(-X));

```

just use

```plaintext
k = inv_logit(X);

```

We’ve gone to great lengths to make our built-ins more numerically stable.

If you define a transformed data variable

```stan
transformed data {
   row_vector[9] exponents = [0, 1, 2, 3, 4, 5, 6, 7, 8];
}

```

you can simplify the body of the function `G` to

```plaintext
return exp(sum(beta_G^exponents));

```

Or you could simplify and just have

```plaintext
return exp(beta_G^[0, 1, 2, 3, 4, 5, 6, 7, 8]);

```

I’d also suggest simplifying the system ode to:

```plaintext
return [G_tracer - k * y[1], G - k * y[2] ]';

```

but I don’t see where `k` is coming from and there is a type error in that `y` needs to be a container if you are going to index it.

---

<div class="post-metadata">

**Author:** ![mbrabec](https://avatars.discourse-cdn.com/v4/letter/m/d6d6ee/32.png) [@mbrabec](https://discourse.mc-stan.org/u/mbrabec)\
**Post date:** [November 2, 2023, 7:47pm UTC](https://discourse.mc-stan.org/t/ode-with-time-varying-parameters-whose-trajectory-is-just-polynomial/33090/3 "2023-11-02T19:47:50Z")

</div>

Thanks a lot!  
I will surely remember not to use named arguments any more …

I fixed the named arguments problem, used inv\_logit, avoided renaming all the elements of arrays in real k and real G definitions and corrected the y to be a vector. Then I get admittedly shorter and more readable piece I list below.

With that, I get further down the line, upt to the error message:  
Error in stanc(file = file, model\_code = model\_code, model\_name = model\_name, :  
0

Semantic error in ‘string’, line 9, column 24 to column 32:

Ill-typed arguments supplied to infix operator \*. Available signatures:  
(int, int) =\> int  
(real, real) =\> real  
(row\_vector, vector) =\> real  
(real, vector) =\> vector  
(vector, real) =\> vector  
(matrix, vector) =\> vector  
(real, row\_vector) =\> row\_vector  
(row\_vector, real) =\> row\_vector  
(row\_vector, matrix) =\> row\_vector  
(real, matrix) =\> matrix  
(vector, row\_vector) =\> matrix  
(matrix, real) =\> matrix  
(matrix, matrix) =\> matrix  
Instead supplied arguments of incompatible type: (real, vector, row\_vector) =\> real, real.

My impression is that this points to some flaw in my construction of the time-varying parameters to be functions, is not it?  
The intention is to solve:

y’\_1(t) = G\_tracer - k(t) \* y\_1(t)  
y’\_2(t) = G(t) - k(t) \* y\_2(t)

where k and G are defined functions.

When I replace the functions k(t) and G(t) by constants in the code,  
i.e. when I replace: return [G\_tracer - k \* y[1], G - k \* y[2] ]‘;  
by: return [G\_tracer - 1 \* y[1], 1 - 1 \* y[2] ]’;  
I get further  
to get error later in calling ode\_rk45:  
Error in stanc(file = file, model\_code = model\_code, model\_name = model\_name, :  
0

Semantic error in ‘string’, line 34, column 30 to column 102:

Ill-typed arguments supplied to function ‘ode\_rk45’. Expected arguments:  
(real, vector) =\> vector, vector, real, real  
Instead supplied arguments of incompatible type:  
(real, vector, vector, vector, real) =\> vector, vector, real, vector, vector, vector, real

P.S.  
I am sorry that I did not really get the return exp(beta\_G\[1\]); idea.  
Somehow, I am missing the t argument there.  
But that is a bit collateral to the main problem, so at the moment I keep the unpleasant but working notation inside the inv\_logit or exp functions.

functions {  
real k(real t, vector beta\_k, row\_vector exponents){  
return inv\_logit(beta\_k[1] + beta\_k[2] \* t + beta\_k[3] \* t^2 + beta\_k[4] \* t^3 + beta\_k[5] \* t^4 + beta\_k[6] \* t^5 +beta\_k[7] \* t^6 +beta\_k[8] \* t^7 + beta\_k[9] \* t^8);  
}  
real G(real t, vector beta\_G){  
return exp(beta\_G[1] + beta\_G[2] \* t + beta\_G[3] \* t^2 + beta\_G[4] \* t^3 + beta\_G[5] \* t^4 + beta\_G[6] \* t^5 + beta\_G[7] \* t^6 + beta\_G[8] \* t^7 + beta\_G[9] \* t^8);  
}  
vector sysode(real t, vector y, vector beta\_k, vector beta\_G, real G\_tracer) {  
return [G\_tracer - k \* y[1], G - k \* y[2] ]';  
}  
}  
data {  
int \< lower =0\> n\_tracer;  
int \< lower =0\> n\_radon;  
vector [n\_tracer] tracer;  
vector [n\_radon] rn;  
real start\_time;  
vector [n\_tracer] times\_tracer;  
}

parameters {  
real\<lower=0\> sigma\_tracer;  
real\<lower=0\> sigma\_rn;  
real\<lower=0\> G\_tracer;  
vector [9] beta\_k;  
vector [9] beta\_G;  
vector [2] y0;  
}  
transformed parameters {  
vector \<lower=0\> [n\_tracer] l\_tracer;  
vector \<lower=0\> [n\_tracer] l\_rn;  
vector \<lower=0\> [n\_radon] mu\_rn;  
array[n\_tracer] vector[2] y = ode\_rk45(sysode, y0, start\_time, times\_tracer, beta\_k, beta\_G, G\_tracer);  
l\_tracer=y[,1];  
l\_rn=y[,2];  
for(j in 2:(n\_radon-1)){  
mu\_rn[j]=sum(l\_rn[((j-2)\*12+2):((j-1)\*12+1)]);  
//mu\_rn[j]=sum(l\_rn[((j-2)\*12+2):((j-1)\*12+1)])/12;  
}  
mu\_rn[1]=mu\_rn[2];  
mu\_rn[n\_radon]=sum(l\_rn[((n\_radon-2)\*12+2):((n\_radon-1)\*12)]);  
//mu\_rn[n\_radon]=sum(l\_rn[((n\_radon-2)\*12+2):((n\_radon-1)\*12)])/12;  
}  
model {  
G\_tracer ~ normal(0, 10);  
beta\_k~normal(0,1);  
beta\_G~normal(0,10);  
sigma\_tracer ~ normal(0, 10);  
sigma\_rn ~ normal(0, 10);  
y0[1]~lognormal(log(50),log(50));  
y0[2]~lognormal(log(500),log(500));

for(i in 1:n\_tracer){  
tracer[i] ~ normal (l\_tracer[i], sigma\_tracer);  
}  
for(j in 1:n\_radon){  
rn[j] ~ normal (mu\_rn[j], sigma\_rn);  
}  
}

* * *

1. 0, 1, 2, 3, 4, 5, 6, 7, 8