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

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.

``````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;
return dy_dt;
}
}
data {
int < lower =0> n_tracer;
int < lower =0> n_radon;
vector [n_tracer] tracer;
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];
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];
}
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);
}
rn[j] ~ normal (mu_rn[j], sigma_rn);
}
}
``````

[edit: escaped Stan code and auto-indented]

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

``````  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

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

just use

``````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

``````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

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

Or you could simplify and just have

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

I’d also suggest simplifying the system ode to:

``````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.

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