Error in _rng function

Modeling
1

Hi
I am trying to fit an ODE system to synthetic data using the following stan code and getting an error for generated quantities block when i fit the model. I am not sure why the error is given although the syntax i believe is correct for normal random number generation. Please help

The error says:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for:

normal_rng(real, int)

Available argument signatures for normal_rng:

normal_rng(real, real)

error in ‘modelc848432fff_myModel_fit’ at line 93, column 56

91:                      
92:                        for (t in 1:T)
93:                          y_ppc[t] = normal_rng(y[t], 1);
                                                           ^
94:                          // y_ppc[t,2] = normal_rng(y[t,2], 1);//this is incorrect

Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model ‘myModel-fit’ due to the above error.

//This is stan code for fitting the model to fake data. The file is called "myModel-fit.stan"
functions {
     real[] ode(real t,
     real[] y,
     real[] params,
     real[] x_r, // x_r:real data, x_i: integer data
     int[] x_i) {

     real dydt[6];
    
     real A;
     real beta;
     real delta;
     real gamma;
     real omega;
     real sigma;
     real kappa;
     real eta;
     delta =params[1]; // death rate of cells infected with virus(this is estimated first)
     A = x_r[1];
     beta = x_r[2];
     gamma = x_r[3];
     omega = x_r[4];
     sigma = x_r[5];
     kappa = x_r[6];
     eta = x_r[7];
     
     dydt[1] =  A - beta * y[5] * y[1] - beta * y[6] * y[1];                 
     dydt[2] = beta * y[1] * y[5] - beta * y[6] * y[2]- delta * y[2];    
     dydt[3] = beta * y[6] * y[1] - beta * y[5] * y[3] ;                 
     dydt[4]=  beta * y[6] * y[2] +  beta * y[5] * y[3] - gamma * y[4];      
     dydt[5] = omega * y[2] + omega * (1-sigma) * y[4] - kappa * y[5];        
     dydt[6] = eta * sigma* y[4] - kappa * y[6];                          
    
  return dydt;
 }
}

data {
    int<lower=1> T; //max time
     real y0[6];   //initial condition of ode
     real t0;      //initial time
     real ts[T];    //time span
     real A;
     real beta;
     real gamma;
     real omega;
     real sigma;
     real kappa;
     real eta;
     real y_hat[T,6];
    // real s; //standatd deviation used in the generated quantities 
}

transformed data {
      real x_r[7];
      int x_i[0];
     
      x_r[1] = A;
      x_r[2] = beta;
      x_r[3] = gamma;
      x_r[4]= omega;
      x_r[5]= sigma;
      x_r[6] = kappa;
      x_r[7] = eta;
}
//the following parameters are to be estimated
parameters{
        real<lower =0> delta;
}
transformed parameters {
          real y[T, 6];
              {
                real params[1];
                
                  params[1] = delta;
                 
                  y = integrate_ode_bdf(ode, y0, t0, ts, params, x_r, x_i);
             }
}
model {
  
  delta ~ normal(0, 10);//prior for delta (this gives delta a half-normal prior as it is declared as lower=0)
  
  for (t in 1:T)
  y_hat[t] ~ normal(y[t], 1);
}
generated quantities {
                     real y_ppc[T, 6];
                     
                       for (t in 1:T)
                         y_ppc[t] = normal_rng(y[t], 1);
                       
}

The data I used is generated from the ODE with some noise as follows

//This is stan code for simulating fake data. The file is called "myModel-sim.stan"

  functions {
     real[] ode(real t,
     real[] y,
     real[] params,
     real[] x_r, // x_r:real data, x_i: integer data
     int[] x_i) {

     real dydt[6];
    
     real A;
     real beta;
     real delta;
     real gamma;
     real omega;
     real sigma;
     real kappa;
     real eta;
     delta =params[1]; // death rate of cells infected with virus(this is estimated first)
     A = x_r[1];
     beta = x_r[2];
     gamma = x_r[3];
     omega = x_r[4];
     sigma = x_r[5];
     kappa = x_r[6];
     eta = x_r[7];
     
     dydt[1] =  A - beta * y[5] * y[1] - beta * y[6] * y[1];                 //Susceptible cells ;
     dydt[2] = beta * y[1] * y[5] - beta * y[6] * y[2]- delta * y[2];      //virus infected cells 
     dydt[3] = beta * y[6] * y[1] - beta * y[5] * y[3] ;                   //TIP infected  cells 
     dydt[4]=  beta * y[6] * y[2] +  beta * y[5] * y[3] - gamma * y[4];        //Co-infected cells
     dydt[5] = omega * y[2] + omega * (1-sigma) * y[4] - kappa * y[5];          //free Virus
     dydt[6] = eta * sigma* y[4] - kappa * y[6];                          //free TIPs
    
  return dydt;
 }

}

data {
    int<lower=1> T; //max time
     real y0[6];   //initial condition of ode
     real t0;      //initial time
     real ts[T];    //time span
     real A;
     real beta;
     real delta;
     real gamma;
     real omega;
     real sigma;
     real kappa;
     real eta;
}

transformed data {
      real params[1];
      real x_r[7];
      int x_i[0];
     
      params[1] = delta;
      x_r[1] = A;
      x_r[2] = beta;
      x_r[3] = gamma;
      x_r[4]= omega;
      x_r[5]= sigma;
      x_r[6] = kappa;
      x_r[7] = eta;
}
model {
  
}
generated quantities {
     real y_hat[T,6];
     y_hat = integrate_ode_bdf(ode, y0, t0, ts, params, x_r, x_i);
// generate fake data and add measurement error
   for (t in 1:T) {
   y_hat[t,1] = y_hat[t,1] + normal_rng(0, 0.1);
   y_hat[t,2] = y_hat[t,2] + normal_rng(0, 0.1);
   y_hat[t,3] = y_hat[t,3] + normal_rng(0, 0.1);
   y_hat[t,4] = y_hat[t,4] + normal_rng(0, 0.1);
   y_hat[t,5] = y_hat[t,5] + normal_rng(0, 0.1);
   y_hat[t,6] = y_hat[t,6] + normal_rng(0, 0.1);
  }
}

And the R code used to generate synthetic data is

sim_data <- list(T = 10, y0 = c(1e8, 0, 0, 0, 1, 1), t0 = 0,
                 ts = 1:10,
                 A= 1.4e6, beta = 3.8e-11, delta =3.3, gamma = 0.14, omega = 1e4, sigma = 0.5, kappa = 3.5, eta = 20000)

#simulate fake data using the system of ODEs
fit1 <- stan("myModel-sim.stan", data = sim_data,
             chains = 1, iter = 1,
             algorithm = "Fixed_param")
#extract data 

y <- extract(fit1, pars = "y_hat")$y_hat[1,,]

The data used to fit the model is

stan_data <- list(T = 10, y0 = c(1e8, 0, 0, 0, 1, 1), t0 = 0,  ts = 1:10,
                 A= 1.4e6, beta = 3.8e-11,  gamma = 0.14, omega = 1e4, 
                 sigma = 0.5, kappa = 3.5, eta = 20000, y_hat = y)
fit2 <- stan("myModel-fit.stan", data = stan_data) #this works fine
2

Per the error message, you have an int type where you need a real type. Also y is a T \times 6 array, and the normal_rng function is not vectorised. Change the generated quantities block to:

generated quantities {
  real y_ppc[T, 6];
  
  for (t in 1:T) {
    for (z in 1:6) {
      y_ppc[t, z] = normal_rng(y[t, z], 1.0);
    }
  }
                       
}

and you should be all good.

3

Thank you very much @hhau. It works

I got another question which is related to the ‘myModel-fit’ code.
I would like to change the model block so that log of y_hat[,5] is normal with mean _log(y[,5])and standard deviation 1. Other vectors in _ y_hat are not modeled.

In Stan, would this be like (after removing the for() loop):
y_hat[,5] ~ lognormal(log(y[,5]), 1);

Thank you

4

Hi again,
When I changed the model and generated quantities block in the above ‘myModel-fit’, I am getting the following error. Could anyone please help. I have already set a lower bound to the parameter, delta.

SAMPLING FOR MODEL 'myModel-fit' NOW (CHAIN 1).
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.00592358, but must be >= 0!  (in 'modela64d712043_myModel_fit' at line 87)

and the last part says

Initialization between (-2, 2) failed after 100 attempts. 
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."

The line 87 is

 y_hat[t] ~ lognormal(log(y[t]), 1);

Here are the code blocks changed.

model {//version 2
  
  delta ~ normal(0, 10);//prior for delta (this gives delta a half-normal prior as it is declared as lower=0)
  
  for (t in 1:T)
  y_hat[t] ~ lognormal(log(y[t]), 1); //version 1
 }
 generated quantities {//version 2:compute predictions from ODE model
       real y_ppc[T, 6];
  
              for (t in 1:T) {
                        for (z in 1:6) { //version 1
                         y_ppc[t, z] = lognormal_rng(log(y[t, z]), 1.0); //version 1
                  }
                }
 }

Thank you

5

Could the problem be that y_hat[2] = -0.00592358?

6

@stijn , I am not sure if the problem is due to y_hat[2] = -0.0059.
Could someone please help me to fix this problem.

Thanks

7

The error message says that you have a negative number in the lognormal statement at line 87. Based on the line numbers in your very first reported error in the original post, I suspect that line 87 is

y_hat[t] ~ lognormal(log(y[t]), 1);

I might be wrong but if I am right you want to check whether all values in y_hat are positive. You can do that in R when you generate the data. You can also make sure that Stan checks whether y_hat is positive by declaring y_hat as positive.

data{
   real<lower=0> y_hat[T,6];
}
8

@stijn Thanks. The stan program now runs okay after changing all values of y_hat positive in R and then declaring y_hat as

data{
   real<lower=0> y_hat[T,6];
}

in Stan program.

1 Like