Row_vector operation issue

Modeling
1

Hi all,

The following stan code, is working. However, when I call it from R, I get the following error regarding to the parameters pAB, uAB, and uASC that are found in the transformed parameters:

[2] " Exception: multiply: Columns of rv (1) and Rows of v (60) must match in size (in ‘string’, line 33, column 2 to column 18)"
[1] “error occurred during calling the sampler; sampling not done”

pAB0, needs to have only one value, hence [1], and rpAB needs to have the dimension of each individual. I based this multiplication on: 4.1 Reductions | Stan Functions Reference
There, I can define a real that can multiply a row_vector and a vector[]. However, a row_vector alone does not work and I need to specify a dimension. Any idea why this error?

functions {
  vector model3(real t, vector y, real pAB, real uAB, real 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 {
  row_vector[1] pAB0;
  row_vector[1] uAB0;
  row_vector[1] uASC0;
  real <lower=0> sigma;
  real <lower=0> sigmapAB;
  real <lower=0> sigmauAB;
  real <lower=0> sigmauASC;
  vector[indivs] rpAB;
  vector[indivs] ruAB;
  vector[indivs] ruASC;
}
transformed parameters {
  real pAB;
  pAB = pAB0*rpAB;
  real uAB;
  uAB = uAB0*ruAB;
  real uASC;
  uASC = uASC0*ruASC;
  array[nobs] vector[2] yhat = ode_bdf(model3, y0, t0, ts, pAB, uAB, uASC);
}
model {
  rpAB ~ lognormal(-2, sigmapAB);
  ruAB ~ lognormal(-1, sigmauAB);
  ruASC ~ lognormal(-0.1, sigmauASC);
  sigmapAB ~ gamma(0.001, 0.001);
  sigmauAB ~ gamma(0.001, 0.001);
  sigmauASC ~ 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); 
  }
}
2

pAB0 is a row vector of length one. rpAB is a (column) vector of length indivs (from the error message we can see that it is of length 60; i.e. indivs = 60).

Stan treats the multiplication of a row vector times a column vector like a matrix multiplication (i.e. a dot product), and thus the number of rows in the column vector must equal the number of columns in the row vector.

You’ve also declared pAB to be a real, and this makes me a bit unsure of what you want.

If you really want pAB to be a vector of length indivs, then you need to declare it as a vector. Then you should declare pAB0 to be a real. To support your intuition here, note that it would also work to leaving pAB0 declared as a 1-dimensional row vector and use pAB0[1] in place of pAB0 at line 33.

If you what you actually want is for pAB to be a real equal to the dot product of rpAB and an indivs-dimensional vector whose elements are all equal to pAB0[1], then you could pass pAB0 as an indivs-dimensional row vector or alternatively you could declare real pAB0 and then do pAB0 * sum(rpAB)

3

Hi,

Thanks for the message. pAB should be a real; however, it has to be transformed by a random effect with indivs dimension and a real parameter. I have tried both options, unfortunately, sum(rpAB) declaring pAB as a real did not give me the results that I was looking for; therefore I declared this parameter as a vector, with the following code, when I get the following error message which I do not really know how to fix it:

Semantic error in ‘string’, line 41, column 31 to column 74:

Ill-typed arguments supplied to function ‘ode_bdf’. Expected arguments:
(real, vector) => vector, vector, real, real[]
Instead supplied arguments of incompatible type:
(real, vector, vector, vector, vector) => vector[], vector, real, real[], vector, vector, vector

functions {
  vector[] model3(real t, vector y, vector pAB, 
  vector uAB, vector uASC) {
    
    vector[rows(y)] dAB_dt = pAB*y[2]-uAB*y[1];
    vector[rows(y)] dASC_dt = -uASC*y[2];
    vector[rows(y)] out[2];
    out = {dAB_dt, dASC_dt};
    
   return (out);
  }
}
data {
  int <lower=1> nobs;
  real t0;
  vector[2] y0;
  real ts[nobs];
  int <lower=1> indivs;
  real <lower=0> antib[nobs];
  //real <lower=1, upper=indivs> subj[nobs];
}
parameters {
  row_vector[1] pAB0;
  row_vector[1] uAB0;
  row_vector[1] uASC0;
  real <lower=0> sigma;
  //real <lower=0> sigmapAB;
  //real <lower=0> sigmauAB;
  //real <lower=0> sigmauASC;
  vector[indivs] rpAB;
  vector[indivs] ruAB;
  vector[indivs] ruASC;
}
transformed parameters {
  vector[indivs] pAB;
  pAB = pAB0[1]*rpAB;
  vector[indivs] uAB;
  uAB = uAB0[1]*ruAB;
  vector[indivs] uASC;
  uASC = uASC0[1]*ruASC;
  array[nobs] vector[2] yhat = ode_bdf(model3, y0, t0, ts, pAB, uAB, uASC);
}
model {
  rpAB ~ lognormal(-2, sqrt(4));
  ruAB ~ lognormal(-1, sqrt(2));
  ruASC ~ lognormal(-0.1, 0.44);
  //sigmapAB ~ gamma(4, 2);
  //sigmauAB ~ gamma(3, 1.5);
  //sigmauASC ~ gamma(1, 2);
  sigma ~ normal(0, 1);
  for (i in 1:nobs) {
    antib[i] ~ lognormal(yhat[i,2], sigma); 
  }
  //antib ~ lognormal(yhat[ : , 2], sigma);
  
}
generated quantities {
  real z_pred[nobs];
  for (j in 1:nobs) {
    z_pred[j] = lognormal_rng(yhat[j,2], sigma); 
  }
}

Edited by @jsocolar for syntax highlighting

4

The error here is from the definition of your ODE function model3, according to the documentation the return type must be vector, but you’re returning an array of vectors (vector[]).

Here’s the section on how to specify the ODE function which will have some more specific tips: 10.2 Ordinary differential equation (ODE) solvers | Stan Functions Reference

1 Like
5

Thanks for the reply. I have noticed that if I would remove the [] from the vector statement in the model3, the code would not work because otherwise I get an error message saying that the left hand side is of type real?

6

It’s more than just removing the [] from the statement, you need to change the values being returned so that it’s a single vector, not an array of vectors.

What is the model code that you’re working with now?

7

All right, thanks! I was using the exact same code as specified above.