Fitting an ODE model; Rejecting initial value

Modeling
,
1

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[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
  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
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_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[1] 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[1] 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[1] 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.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 500 [  0%]  (Warmup)
2 Likes
2

You define alpha to be non-negative, but that allows values greater than 1. As a result, alpha[6] 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.

1 Like
3

That makes sense––thanks!