ODE with time-varying parameters whose trajectory is just polynomial

mbrabec · October 18, 2023, 3:38pm

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

Bob_Carpenter · November 1, 2023, 9:20pm

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.

mbrabec · November 2, 2023, 7:47pm

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

0, 1, 2, 3, 4, 5, 6, 7, 8 ↩︎

Topic		Replies	Views
Dealing with unknown data with rstan/ODE Modeling rstan , ode	2	431	November 10, 2020
ODEs with time dependent parameters Modeling	7	2962	November 9, 2017
ODE vector function Modeling	2	766	March 14, 2022
Updating ode to new syntax. Exception: ode_rk45: ode paramers and data is inf Modeling rstan , fitting-issues , ode , cmdstanr	1	50	July 17, 2024
Using ODE solver with varying upper and lower bounds Modeling rstan , specification , ode	6	883	April 11, 2020

ODE with time-varying parameters whose trajectory is just polynomial

Related topics