 # Fitting an ODE model; Rejecting initial value

Hey guys, I’m trying to fit a model in Stan than involves an ODE model. I keep getting a “Rejecting initial value” error that I’m not sure how to fix. I’m very new to this, so any help would be greatly appreciated.

For context, I’m working in rstan.

``````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; // System state coded as an array, such that Z = (P,H)
real H = Z;

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

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; // Initial conditions for ODE
int<lower=0>y[N,2]; // Define y as real and non-negative
}

parameters{
real<lower=0>alpha; // Make all items in alpha non-neg
real<lower=0>Z_init; // 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
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
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_Data//TypeII_Fitting_Output2.R",
seed = 123,
show_messages = TRUE,
verbose = TRUE)
``````

And this is what the output 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: beta_lpdf: Random variable is 2.07072, but must be less than or equal to 1  (in 'model6c63eadfb23_Stan_Model_TypeII' at line 45)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: beta_lpdf: Random variable is 2.03603, but must be less than or equal to 1  (in 'model6c63eadfb23_Stan_Model_TypeII' at line 45)

Chain 1: Rejecting initial value:
Chain 1:   Error evaluating the log probability at the initial value.
Chain 1: Exception: beta_lpdf: Random variable is 5.10664, but must be less than or equal to 1  (in 'model6c63eadfb23_Stan_Model_TypeII' at line 49)

Chain 1:
Chain 1: Gradient evaluation took 0.010238 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 102.38 seconds.
You define `alpha` to be non-negative, but that allows values greater than 1. As a result, `alpha` can be outside (0, 1), which makes its probability under the Beta prior be zero. You can define each parameter separately and then wrap them up in a vector/array `theta` or something in order to pass to the ODE solver.