ODE with parameters as vectors

Modeling
1

Hi all,

I was wondering to use a function like mentioned below, but instead of using theta as a real[] to use instead vector theta? I would like to transform later on in the code these parameters of theta, and in order to transform them they need to be a vector type, but the problem is that I have already define this a real. Would that be possible?
Do I also need then to change the type real model3 to vector model3 as well?

functions {
  real [] model3(real t, real[] y, real[] theta,
  real[] x_r, int[] x_i) {
    
    real AB= y[1];
    real ASC = y[2];
    
    real pAB = theta[1];
    real uAB = theta[2]; 
    real uASC = theta[3];
    
    
    real dAB_dt = pAB*ASC - uAB*AB;
    real dASC_dt = -uASC*ASC;
    
   return {dAB_dt, dASC_dt};
  }
}
2

use to_vector, see 8 Mixed Operations | Stan Functions Reference

1 Like
3

Thanks!
However, I have tried this but since I defined in the function these parameters as reals it did not work; I have changed the code with the following code, with the following error message:
Error in stanc(filename, allow_undefined = TRUE) : 0

Semantic error in ‘string’, line 7, column 4 to column 32:

Ill-typed arguments supplied to assignment operator =: lhs has type real and rhs has type vector

functions {
  vector model3(real t, vector y, vector pAB, 
  real uAB, vector uASC) {
    
    vector[2] dydt;
    
    dydt[1] = pAB*y[2]-uAB*y[1];
    dydt[2] = -uASC*y[2];
    
   return dydt;
  }
}
data {
  int <lower=1> nobs;
  real t0;
  vector[2] y0;
  array[nobs] real ts;
  int <lower=1> indivs;
  array[nobs] real antib;
  real <lower=1, upper=indivs> subj[nobs];
}
parameters {
  //real <lower=0> AB0;
  real <lower=0> pAB0;
  real <lower=0> uAB;
  real <lower=0> uASC0;
  //real <lower=0> sigmaAB0;
  real <lower=0> sigma;
  vector[indivs] rpAB;
  vector[indivs] ruASC;
}
transformed parameters {
  vector[indivs] pAB;
  pAB = pAB0*rpAB;
  vector[indivs] uASC;
  uASC = uASC0*ruASC;
  array[nobs] vector[2] yhat = ode_bdf_tol(model3, y0, t0, ts, 1e-8, 1e-8, 1000, pAB, uAB, uASC);
}
model {
  rpAB ~ lognormal(-2, sqrt(4));
  uAB ~ lognormal(-1, sqrt(2));
  ruASC ~ lognormal(-0.1, 0.44);
  //AB0 ~ lognormal(-2, 4);
  //sigmaAB0 ~ gamma(0.001, 0.001);
  sigma ~ normal(0, 1);
  
  //for (i in 1:nobs) {
    //antib[i] ~ lognormal(yhat[i,1], sigma); 
  //}
  antib ~ lognormal(log(yhat[ : , 2]), sigma);
}
generated quantities {
  array[nobs] real z_pred;
  for (j in 1:nobs) {
    z_pred[j] = lognormal_rng(log(yhat[j,2]), sigma); 
  }
}

Do you know why I could not use this vector type in the function denoted as pAB and uASC?

4

If a function takes a real [] argument and you for some reason have the parameters as vector (or others like row_vector or a 1D matrix) you should use to_array_1d to convert them to the right type.

This will not change the type of the function itself, which will remain real []; however, note that for ODEs this is the old interface, the new interface has output type vector and takes parameters as separate real arguments (and I guess can also take a packed real array just like the old one).

1 Like
5

The error arises because at line 7 dydt[1] has type real. (dydt is a vector; its elements are reals). pAB*y[2]-uAB*y[1]; has type vector, because pAB is a vector and a vector times a scalar is also a vector.

Thus, we see the error

Semantic error in ‘string’, line 7, column 4 to column 32:

Ill-typed arguments supplied to assignment operator =: lhs has type real and rhs has type vector

It’s telling us that the left hand side is a real, and the right hand side is a vector.

1 Like
6

Where do you mean that I should use the argument to_array_1d?
I have tried also the following in the transformed parameters, defining in the function everything as real (for the parameters) instead of vectors:

to_vector(pAB) = pAB0*rpAB

However, this did not worked either

7

Thanks! However, when I change uAB to a vector type, and I also do the same transformation for uAB in the transformed parameters I do still get the same error message, please find the code:

functions {
  vector model3(real t, vector y, vector pAB, 
  vector uAB, vector uASC) {
    
    vector[2] dydt;
    
    dydt[1] = pAB*y[2]-uAB*y[1];
    dydt[2] = -uASC*y[2];
    
   return dydt;
  }
}
data {
  int <lower=1> nobs;
  real t0;
  vector[2] y0;
  array[nobs] real ts;
  int <lower=1> indivs;
  array[nobs] real antib;
  real <lower=1, upper=indivs> subj[nobs];
}
parameters {
  //real <lower=0> AB0;
  real <lower=0> pAB0;
  real <lower=0> uAB0;
  real <lower=0> uASC0;
  //real <lower=0> sigmaAB0;
  real <lower=0> sigma;
  vector[indivs] rpAB;
  vector[indivs] ruAB;
  vector[indivs] ruASC;
}
transformed parameters {
  vector[indivs] pAB;
  pAB = pAB0*rpAB;
  vector[indivs] uAB;
  uAB = uAB0*ruAB;
  vector[indivs] uASC;
  uASC = uASC0*ruASC;
  array[nobs] vector[2] yhat = ode_bdf_tol(model3, y0, t0, ts, 1e-8, 1e-8, 1000, pAB, uAB, uASC);
}
model {
  rpAB ~ lognormal(-2, sqrt(4));
  ruAB ~ lognormal(-1, sqrt(2));
  ruASC ~ lognormal(-0.1, 0.44);
  //AB0 ~ lognormal(-2, 4);
  //sigmaAB0 ~ gamma(0.001, 0.001);
  sigma ~ normal(0, 1);
  
  //for (i in 1:nobs) {
    //antib[i] ~ lognormal(yhat[i,1], sigma); 
  //}
  antib ~ lognormal(log(yhat[ : , 2]), sigma);
}
generated quantities {
  array[nobs] real z_pred;
  for (j in 1:nobs) {
    z_pred[j] = lognormal_rng(log(yhat[j,2]), sigma); 
  }
}
´´´
8

The function arguments tell us that pAB is a vector, and uAB is a vector. At line 7, pAB*y[2]-uAB*y[1] is therefore a vector. But at line 7 dydt[1] is not a vector; it is a real.

1 Like
9

Thanks, this helped

1 Like