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
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)"
[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.