Can't concatenate two vectors to fix "Rejecting initial value" error

Hey all! I had an issue yesterday, wherein I was getting a “Rejecting initial value” error that wasn’t letting my model run. I have log-normally and beta distributed parameters, which I put a lower bound of 0 on, and no upper bound. I didn’t include an upper bound on the beta distributed parameters because I couldn’t figure out how to put different bounds on different parameters. Unfortunately, I was told that this was where my “Rejecting initial values” error was coming from. So, now I’m trying to split my parameters into two vectors, which I can combine into one to pass through the ODE solver, and I have no idea what I’m doing.

Any help would be greatly appreciated!

Here’s my original script with the “Rejecting initial values” error:

library(rstan)
library(gdata)
library(bayesplot)
library(tidyverse)
library(brms)

write("// Stan Model;
functions{
  real[] dZ_dt(
    real t, // Time
    real[] Z, // System state {Parasite, Host}
    real[] alpha, // Parameters 
    real[] x_r, // Real data; below is integer data
    int[] x_i){
    real P = Z[1]; // System state coded as an array, such that Z = (P,H)
    real H = Z[2];

    real r = alpha[1]; // Parameters of the system, in order they appear
    real O = alpha[2];
    real h = alpha[3];
    real b = alpha[4];
    real c = alpha[5];
    real u = alpha[6];

    real dP_dt = P*r - H*(O*P/1 + O*P*h); // Mechanistic model
    real dH_dt = b + H*(c*(O*P/1 + O*P*h)-u);
    return{dP_dt,dH_dt}; // Return the system state
  }
}

data{
  int<lower=0>N; // Define N as non-negative integer
  real ts[N]; // Assigns time points to N
  int y_init[2]; // Initial conditions for ODE
  int<lower=0>y[N,2]; // Define y as real and non-negative
}

parameters{
  real<lower=0>alpha[6]; // Make all items in alpha non-neg
  real<lower=0>Z_init[2]; // Initial population size non-neg
}

transformed parameters{
  // Stiff solver for ODE
  real Z[N,2]
  = integrate_ode_bdf(dZ_dt,Z_init,0,ts,alpha,rep_array(0.0,0),rep_array(0,0),1e-6,1e-5,2000);
}

model{ // Ignore the means/variances for now...
  alpha[{1}]~lognormal(log(2.7),1); // r
  alpha[{2}]~beta(2,2); // O; bounded between 0 and 1
  alpha[{3}]~lognormal(log(10),1); // h
  alpha[{4}]~lognormal(log(18),1); // b
  alpha[{5}]~lognormal(log(0.02),1); // c
  alpha[{6}]~beta(2,2); // u; bounded between 0 and 1
  Z_init~lognormal(log(10),1);
  for (k in 1:2){
    y_init[k]~poisson(Z_init[k]);
    y[ ,k]~poisson(Z[ ,k]);
  }
}",
"Stan_Model_TypeII.stan")

stanc("Stan_Model_TypeII.stan") # To check that we wrote a file (I think we did?)

Stan_Model_TypeII <- stan_model("Stan_Model_TypeII.stan")

# Squeezing the data into a form that Stan gets
Stoch_Data_TypeII = read_csv('/Users/mjarviscross/Desktop/GitHub/FcnalRespProj/data/Type II/FcnalRespProj_Data/Stoch_Data_TypeII.csv')
N <- length(Stoch_Data_TypeII$t)-1 # N is 1952 which makes sense bc length of DF is 1953
ts <- 1:N
y_init <- c(1,18) # Initial states, P = 1; H = 18
y <- as.matrix(Stoch_Data_TypeII[2:(N+1),2:3])
y <- cbind(y[ ,2],y[ ,1]); # This worked, sick; where y[,1] is H, and y[,2] is P
Stan_StochData_TypeII <- list(N=N,ts=ts,y_init=y_init,y=y)

# // Make non-neg real<lower=0>alpha[{1,3,4,5}];
# // Proportions bounded between 0 and 1 real<lower=0,upper=1>alpha[{2,6}]

# Fitting the data to the model
Stan_Model_TypeII <- stan_model("Stan_Model_TypeII.stan")
fit <- sampling(Stan_Model_TypeII,
                data = Stan_StochData_TypeII,
                warmup = 250,
                iter = 500,
                chains = 1,
                cores = 1,
                thin = 1,
                check_data = TRUE,
                diagnostic_file = "/Users/mjarviscross/Desktop/GitHub/FcnalRespProj/data/Type II/FcnalRespProj/TypeII_Fitting_Output2.R",
                seed = 123,
                show_messages = TRUE,
                verbose = TRUE)

And here’s my new script, wherein I can’t figure out how to concatenate two sets of parameters (just the model part):

functions{
  real[] dZ_dt(
    real t, // Time
    real[] Z, // System state {Parasite, Host}
    real[] alpha, // Parms (r,h,b,c)
    real[] alpha2, // Proportional parms (O and u)
    real[] x_r, // Real data; below is integer data
    int[] x_i){
    real P = Z[1]; // System state coded as an array, such that Z = (P,H)
    real H = Z[2];

    real r = alpha[1]; // Parameters of the system, in order they appear
    real O = alpha2[1];
    real h = alpha[2];
    real b = alpha[3];
    real c = alpha[4];
    real u = alpha2[2];

    real dP_dt = P*r - H*(O*P/1 + O*P*h); // Mechanistic model
    real dH_dt = b + H*(c*(O*P/1 + O*P*h)-u);
    return{dP_dt,dH_dt}; // Return the system state
  }
}

data{
  int<lower=0>N; // Define N as non-negative integer
  real ts[N]; // Assigns time points to N
  int y_init[2]; // Initial conditions for ODE
  int<lower=0>y[N,2]; // Define y as real and non-negative
}

parameters{
  real<lower=0>alpha[4];
  real<lower=0,upper=1>alpha2[2]; // Make all items in alpha non-neg
  real<lower=0>Z_init[2]; // Initial population size non-neg
}

transformed parameters{
  // Stiff solver for ODE
  vector[N] alpha3;
  alpha3 = append_row(alpha,alpha2);
  real Z[N,2]
  = integrate_ode_bdf(dZ_dt,Z_init,0,ts,alpha3,rep_array(0.0,0),rep_array(0,0),1e-6,1e-5,2000);
}

model{ // Ignore the means/variances for now...
  alpha[{1}]~lognormal(log(2.7),1); // r
  alpha2[{2}]~beta(2,2); // O; bounded between 0 and 1
  alpha[{3}]~lognormal(log(10),1); // h
  alpha[{4}]~lognormal(log(18),1); // b
  alpha[{5}]~lognormal(log(0.02),1); // c
  alpha2[{2}]~beta(2,2); // u; bounded between 0 and 1
  Z_init~lognormal(log(10),1);
  for (k in 1:2){
    y_init[k]~poisson(Z_init[k]);
    y[ ,k]~poisson(Z[ ,k]);
  }
}

I gather that now, the issue is that alpha and alpha2 are both real objects, and append_row can’t take two real arguments (only (vector x, vector y), (vector x, real y), or (real x, vector y)), but I’m not sure how to make any of these scenarios work…

I tried to assign parameters individually, but can’t get that to work either, and am getting an error in the ODE solver, due to the way I’ve included the parameters.

functions{
  real[] dZ_dt(
    real t, // Time
    real[] Z, // System state {Parasite, Host}
    real[] x_r, // Real data; below is integer data
    int[] x_i){
    real P = Z[1]; // System state coded as an array, such that Z = (P,H)
    real H = Z[2];

    real r;  // Parameters of the system, in order they appear
    real O;
    real h;
    real b;
    real c;
    real u;

    real dP_dt = P*r - H*(O*P/1 + O*P*h); // Mechanistic model
    real dH_dt = b + H*(c*(O*P/1 + O*P*h)-u);
    return{dP_dt,dH_dt}; // Return the system state
  }
}

data{
  int<lower=0>N; // Define N as non-negative integer
  real ts[N]; // Assigns time points to N
  int y_init[2]; // Initial conditions for ODE
  int<lower=0>y[N,2]; // Define y as real and non-negative
}

parameters{
  real<lower=0>r;
  real<lower=0,upper=1>O;
  real<lower=0>h;
  real<lower=0>b;
  real<lower=0>c;
  real<lower=0,upper=1>u;
  real<lower=0>Z_init[2]; // Initial population size non-neg
}

transformed parameters{
  // Stiff solver for ODE
  real Z[N,2]
  = integrate_ode_bdf(dZ_dt,Z_init,0,ts,{r, O, h, b, c, u},rep_array(0.0,0),rep_array(0,0),1e-6,1e-5,2000);
}

model{ // Ignore the means/variances for now...
  r~lognormal(log(2.7),1); // r
  O~beta(2,2); // O; bounded between 0 and 1
  h~lognormal(log(10),1); // h
  b~lognormal(log(18),1); // b
  c~lognormal(log(0.02),1); // c
  u~beta(2,2); // u; bounded between 0 and 1
  Z_init~lognormal(log(10),1);
  for (k in 1:2){
    y_init[k]~poisson(Z_init[k]);
    y[ ,k]~poisson(Z[ ,k]);
  }
}

Thank you so much!

The problem with second program is that you place the concatenated vector in generated quantities rather than transformed parameters. Generated quantities play no role in the model, the goal of the block is to have quantities derived from the model (such as predictions or log-likelihoods) ready for post processing.

Since you only have six parameters, I recommend you actually name them all and assign priors individually. This will also help with readability.

You can then just do

   integrate_ode_bdf(..., 0, ts, {r, o, h, b, c, u}, ...);

assuming you use the same names as in the dZ_dtblock.

Also, you’re still placing a Beta prior on a parameter with no upper bound.

Oops, I just accidentally reversed alpha and alpha2.

I moved the appending lines into the trans parms block, and I’m still getting an error related to the arguments in append_row.

Also, I haven’t listed the parameters separately because I’m not sure how to do that. I’m very new to Stan.

I tried this, but can’t get it to work either… I feel like I’m missing something…

functions{
  real[] dZ_dt(
    real t, // Time
    real[] Z, // System state {Parasite, Host}
    real[] x_r, // Real data; below is integer data
    int[] x_i){
    real P = Z[1]; // System state coded as an array, such that Z = (P,H)
    real H = Z[2];

    real r;  // Parameters of the system, in order they appear
    real O;
    real h;
    real b;
    real c;
    real u;

    real dP_dt = P*r - H*(O*P/1 + O*P*h); // Mechanistic model
    real dH_dt = b + H*(c*(O*P/1 + O*P*h)-u);
    return{dP_dt,dH_dt}; // Return the system state
  }
}

data{
  int<lower=0>N; // Define N as non-negative integer
  real ts[N]; // Assigns time points to N
  int y_init[2]; // Initial conditions for ODE
  int<lower=0>y[N,2]; // Define y as real and non-negative
}

parameters{
  real<lower=0>r;
  real<lower=0,upper=1>O;
  real<lower=0>h;
  real<lower=0>b;
  real<lower=0>c;
  real<lower=0,upper=1>u;
  real<lower=0>Z_init[2]; // Initial population size non-neg
}

transformed parameters{
  // Stiff solver for ODE
  real Z[N,2]
  = integrate_ode_bdf(dZ_dt,Z_init,0,ts,{r, O, h, b, c, u},rep_array(0.0,0),rep_array(0,0),1e-6,1e-5,2000);
}

model{ // Ignore the means/variances for now...
  r~lognormal(log(2.7),1); // r
  O~beta(2,2); // O; bounded between 0 and 1
  h~lognormal(log(10),1); // h
  b~lognormal(log(18),1); // b
  c~lognormal(log(0.02),1); // c
  u~beta(2,2); // u; bounded between 0 and 1
  Z_init~lognormal(log(10),1);
  for (k in 1:2){
    y_init[k]~poisson(Z_init[k]);
    y[ ,k]~poisson(Z[ ,k]);
  }
}
 

Sorry, I shouldn’t have been so curt. Welcome to the community.

Does this still throw the same initialisation error?

@maxbiostat No worries at at all; thanks so much for all your help!

I can’t run the model to look for the initialisation error because there’s still a syntax error. I added the version of the script wherein I tried to assign the parameters individually to the original post, but in the ODE solver, I’m getting an error that reads:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Wrong signature for function integrate_ode_bdf; first argument must be the name of a function with signature (real, real[ ], real[ ], data real[ ], data int[ ]): real[ ].
 error in 'model1e841bc9b83_Updated_Fitting_Script_MJC3' at line 43, column 105
  -------------------------------------------------
    41:   // Stiff solver for ODE
    42:   real Z[N,2]
    43:   = integrate_ode_bdf(dZ_dt,Z_init,0,ts,{r, O, h, b, c, u},rep_array(0.0,0),rep_array(0,0),1e-6,1e-5,2000);
                                                                                                                ^
    44: }
  -------------------------------------------------
1 Like

You were missing the “parameters” argument, theta.

functions{
  real[] dZ_dt(
    real t, // Time
    real[] Z, // System state {Parasite, Host}
    real [] theta, // parameters
    real[] x_r, // Real data; below is integer data
    int[] x_i){
      
    real P = Z[1]; // System state coded as an array, such that Z = (P,H)
    real H = Z[2];

    real r = theta[1];  // Parameters of the system, in order they appear
    real O = theta[2];
    real h = theta[3];
    real b = theta[3];
    real c = theta[3];
    real u = theta[3];

    real dP_dt = P*r - H*(O*P/1 + O*P*h); // Mechanistic model
    real dH_dt = b + H*(c*(O*P/1 + O*P*h)-u);
    return({dP_dt,dH_dt}); // Return the system state
  }
}

data{
  int<lower=0>N; // Define N as non-negative integer
  real ts[N]; // Assigns time points to N
  int y_init[2]; // Initial conditions for ODE
  int<lower=0>y[N,2]; // Define y as real and non-negative
}

parameters{
  real<lower=0>r;
  real<lower=0,upper=1>O;
  real<lower=0>h;
  real<lower=0>b;
  real<lower=0>c;
  real<lower=0,upper=1>u;
  real<lower=0>Z_init[2]; // Initial population size non-neg
}

transformed parameters{
  // Stiff solver for ODE
  real Z[N,2]
  = integrate_ode_bdf(dZ_dt, Z_init, 0.0, ts, {r, O, h, b, c, u}, rep_array(0.0, 2), rep_array(0, 2), 1e-6, 1e-5, 2000);
}

model{ // Ignore the means/variances for now...
  r~lognormal(log(2.7),1); // r
  O~beta(2,2); // O; bounded between 0 and 1
  h~lognormal(log(10),1); // h
  b~lognormal(log(18),1); // b
  c~lognormal(log(0.02),1); // c
  u~beta(2,2); // u; bounded between 0 and 1
  Z_init~lognormal(log(10),1);
  for (k in 1:2){
    y_init[k]~poisson(Z_init[k]);
    y[ ,k]~poisson(Z[ ,k]);
  }
}

I also fixed the rep_array and the return instruction in dZ_dt.

@maxbiostat Thank so much! That makes sense. That fixed the syntactic error, but now I’m getting an error when I run the model, that reads:

CHECKING DATA AND PREPROCESSING FOR MODEL 'Stan_Model_TypeII' NOW.

COMPILING MODEL 'Stan_Model_TypeII' NOW.

STARTING SAMPLER FOR MODEL 'Stan_Model_TypeII' NOW.

SAMPLING FOR MODEL 'Stan_Model_TypeII' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Rate parameter[5] is -2.11967e-07, but must be >= 0!  (in 'model5e34e184223_Stan_Model_TypeII' at line 58)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Rate parameter[6] is -3.47625e-05, but must be >= 0!  (in 'model5e34e184223_Stan_Model_TypeII' at line 58)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Rate parameter[18] is -1.53267e-07, but must be >= 0!  (in 'model5e34e184223_Stan_Model_TypeII' at line 58)

Chain 1: 
Chain 1: Gradient evaluation took 0.010365 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 103.65 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
Chain 1: Exception: CVode failed with error flag -1  (in 'model5e34e184223_Stan_Model_TypeII' at line 44)

That seems to stem from numerical errors that are beyond my knowledge. @bbbales2 @wds15 Halp?

1 Like

It’s entirely possible for the ODE solver to go negative (even if mathematically the eqs shouldn’t).

Given the tolerances:

1e-6, 1e-5

Seeing numbers below zero on the order of the tolerances kinda makes sense in that regard.

You can either write your ODE on the log scale (\frac{\partial \log y}{\partial t} = \frac{1}{y} \frac{\partial y}{\partial t}), or you might get away with just adding a small number to the output of your ODE to make sure things are positive.

But you could check that this is a tolerance thing by making the tolerances smaller and seeing if the errors about things being below zero stay the same.

3 Likes

Hi @bbbales2, thanks for your input! I’m not really sure how to asses these values and change them such that the ODE solver doesn’t go negative. Do you have any suggestions?

Just make them smaller numbers:

real Z[N,2]
  = integrate_ode_bdf(dZ_dt, Z_init, 0.0, ts, {r, O, h, b, c, u}, rep_array(0.0, 2), rep_array(0, 2), 1e-10, 1e-10, 1e4);

For instance (I also made the max number of steps larger) and then see if these numbers get smaller:

Chain 1: Exception: poisson_lpmf: Rate parameter[5] is -2.11967e-07, but must be >= 0!  (in 'model5e34e184223_Stan_Model_TypeII' at line 58)

@bbbales2 I see––thank you so much. I made the changes you suggested and am getting the same errors. The values are getting smaller, but on different parameters? I’m very confused:

CHECKING DATA AND PREPROCESSING FOR MODEL 'Stan_Model_TypeII' NOW.

COMPILING MODEL 'Stan_Model_TypeII' NOW.

STARTING SAMPLER FOR MODEL 'Stan_Model_TypeII' NOW.

SAMPLING FOR MODEL 'Stan_Model_TypeII' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Rate parameter[8] is -9.02654e-19, but must be >= 0!  (in 'model699381bbbe8_Stan_Model_TypeII' at line 58)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Rate parameter[44] is -1.61356e-18, but must be >= 0!  (in 'model699381bbbe8_Stan_Model_TypeII' at line 58)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: poisson_lpmf: Rate parameter[35] is -1.62491e-12, but must be >= 0!  (in 'model699381bbbe8_Stan_Model_TypeII' at line 58)

Chain 1: 
Chain 1: Gradient evaluation took 0.018809 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 188.09 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)

@bbbales2 I made the values even smaller, and am now getting a really strange output?

CHECKING DATA AND PREPROCESSING FOR MODEL 'Stan_Model_TypeII' NOW.

COMPILING MODEL 'Stan_Model_TypeII' NOW.

STARTING SAMPLER FOR MODEL 'Stan_Model_TypeII' NOW.

SAMPLING FOR MODEL 'Stan_Model_TypeII' NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: CVode failed with error flag -2  (in 'model77d7e3f0a6b_Stan_Model_TypeII' at line 44)

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                                               
[2] "  Exception: CVode failed with error flag -2  (in 'model77d7e3f0a6b_Stan_Model_TypeII' at line 44)"                                                                             
[3] "In addition: Warning message:"                                                                                                                                                  
[4] "In writable_sample_file(diagnostic_file) :"                                                                                                                                     
[5] "  \"/Users/mjarviscross/Desktop/GitHub/FcnalRespProj/data/Type II/FcnalRespProj\" is not writable; use \"/var/folders/l_/xb7lfbls12nb0dwxdrww8gn40000gn/T//Rtmp0blNqe\" instead"
1 Like

If you make the tolerances too small the solver will fail to solve to that tolerance. I assume that’s the second error.

The smaller things in different places makes sense. Lowering the tolerances means you’ll solve the eqs more exactly and be less likely to get negative numbers. Presumably -1.62491e-12 is within a couple multiples of the tolerance of zero? Think of it like that – according to the solver everything is fine to within tolerance.

So you’re probably fine to run with a low tolerance and just add a small amount to the output.

If you add 1e-10 to the output everywhere so everything stays positive, then you’re biasing your solutions up, but presumably that won’t matter in terms of the stats. That make sense?

This makes sense to me––thanks!

I’m confused by this bit––would you mind explaining a bit more?

I’m so sorry I’m being difficult––I really really appreciate your help!

1 Like

Sure thing.

Since your likelihood is something like:

y ~ poisson(z);

z is the output of an ODE solver. It turns out for numerical reasons, the ODE solver is giving you z values that are sometimes negative, and that messes up the likelihood evaluation.

If z > -1e-10, then:

y ~ poisson(z + 1e-10);

will not have similar problems, and presumably shifting the mean parameter of the poisson by 1e-10 isn’t gonna mess up the statistics.

That make sense? You can also rewrite your equations to solve them on the log scale I think, but this is an easier solution if it works.

That makes sense to me. I just tried it, and I got a really similar error, now saying that some values were too high given their bounds; for example, that parameters that have an upper bound of 1 are being estimated as 2 or 5.