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.