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!