"Rejecting initial value" and "Error evaluating the log probability at the initial value" when trying to fit with an ODE

Modeling
, ,
1

Hi all!

I’ve been trying to get this fitting model to work for a while, and I keep getting an error that reads:

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 'modelbaf75d77eda_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 'modelbaf75d77eda_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 'modelbaf75d77eda_Stan_Model_TypeII' at line 58)

etc.

I’ve tried to troubleshoot this problem before, and have gotten some great advice from the community. However, nothing seems to be working.

Here’s the actual script (I’m using rstan):

write("// Stan Model;
functions{
  real[] dZ_dt(
    real t, // Time
    real[] Z, // System state {Parasite, Host}
    real[] theta, // Parms
    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[4];
    real c = theta[5];
    real u = theta[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>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,2),rep_array(0,2),1e-10,1e-10,2e4);
}

model{
  r~normal(2.5,1); // r
  O~beta(2,2); // O; bounded between 0 and 1
  h~normal(0.06,1); // h
  b~normal(35,1); // b
  c~normal(0.2,1); // c
  u~beta(2,2); // u; bounded between 0 and 1
  Z_init~uniform(0,400);
  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

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 # df is 2923; N is 2922
ts <- 1:N
y_init <- c(10,350) # Initial states, P = 10; H = 350
y <- as.matrix(Stoch_Data_TypeII[2:(N+1),3:4])
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 = 124,
                show_messages = TRUE,
                verbose = TRUE)

Previously, my priors that are normally distributed were log-normally distributed. I changed them to normal to discourage the sampler from getting too close to 0.

I’ve also gone back to my data, and re-simulated it to keep initial values and population abundances away from 0.

I have no idea what I’m doing wrong, and would really appreciate any advice people are willing to give me.

Thank you so much!

P.S. My priors look super tight because they are; I’m using my actual parameter values as means in hopes of forcing the model to estimate my parameters correctly.

1 Like
2

The error says that the rate parameter in y[ ,k]~poisson(Z[ ,k]); goes below zero. At late times some components of Z fall to approximately zero but the numerical integrator overshoots a little bit and Z becomes a tiny negative number. The easiest fix is to adjust Z upwards a little.

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,2),rep_array(0,2),1e-10,1e-10,2e4);
  for (i in 1:N) {
    Z[i,1] = Z[i,1] + 1e-10;
    Z[i,2] = Z[i,2] + 1e-10;
  }
}
3 Likes
3

Hi @nhuurre! Thanks for the advice! I added the lines you suggested and increased the added value to 1e-6, and it’s not throwing the error anymore! Figures crossed for the sampling… Thank you so much!! :)

1 Like