Hello - I’m trying to fit a fish passage model. After some exposure to BUGS years ago this is my first jump into STAN. I am having some issues and I’m hoping the community could provide some help.
I am fitting a fish passage model based on dx.doi.org/10.1139/cjfas-2014-0554 (but just for a single year at this point). I am fitting a negative binomial model and it runs (not all parameters look very good though) but the issue I’m having is I can figure out how to include a 'normalizing constant (\psi_sum). In my STAN code I have a ‘25’ in place of where that parameter should go but when I sub the \psi_sum in for the 25 I get the error:
“Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: neg_binomial_2_log_lpmf: Log location parameter[1] is nan, but must be finite! (in ‘modela7fc42ebcdfc_naive_single3’ at line 41)”.
- Adding inits doesn’t seem to solve the issue.
Any help would be appreciated. Thank you.
First, I simulate some data:
####### SIMULATE DATA IN R
r <- 10000 # total run
p <- 250 # peak escapement DAY
sigma <- 10 # error
d <- 200:300 #days
set.seed(111)
psi <- sum(exp(-((d-p)^2)/(2*sigma^2)))
run <- r*(exp(-((d-p)^2)/(2*sigma^2)))/psi
run_error <- rlnorm(length(run),log(run),0.2)
single_sim <- data.frame(julian=d,Count=round(run_error,0))
ggplot(data=single_sim,aes(x=julian,y=Count))+ geom_segment(aes(x = julian,xend=julian,y=0,yend=Count),alpha=0.75)+theme_bw()
#########################################################
# Set up sampling conditions
##################################
##### Set STAN conditions
warmups <- 1000
total_iterations <- 2000
max_treedepth <- 10
n_chains <- 1
n_cores <- 4
adapt_delta <- 0.95
### Format data for stan model
sim_dat <- with(single_sim,list(y=Count,x=julian,N=nrow(sim_run)))
### Run stan model
fit <- stan(file = 'naive_single3.stan', data = sim_dat,
chains = n_chains,
warmup = warmups,
iter = total_iterations,
cores = n_cores,
refresh = 250,
init = list(list(log_r=5,p=230,sigma_arrival=5,reciprocal_phi=.1,r=exp(7),phi=.1,c_hat=rep(10,101),mu=rep(10,101,y_rep=rep(10,101),phi2=rep(0.01,101)))))
library(posterior)
########## STAN MODEL CODE ###########
data {
int N; // Number of observations (101)
int y[N]; // Vector of observations
vector[N] x;// Vector of DOY
}
// The parameters we are going to estimate in our model
parameters {
real<lower=0> p;
real<lower=0.01> log_r;
real<lower=0> sigma_arrival;
real<lower=0> reciprocal_phi;
}
transformed parameters{
real r = exp(log_r);
real phi;
// vector[N] phi2;
// real<lower=0> phi_sum;
vector[N] c_hat;
phi = 1. / reciprocal_phi;
for (i in 1:N){
c_hat[i] = r*exp(-(pow((x[i]-p),2)/(2*pow(sigma_arrival,2))))/(25);
}
//for (i in 1:N) {
// phi2[i] = exp(-(pow((x[i]-p),2)/(2*pow(sigma_arrival,2))));
// }
//phi_sum = sum(phi2);
}
model {
//Priors
reciprocal_phi ~ cauchy(0., 3);
log_r ~ normal(log(10000),5);
y ~ neg_binomial_2_log(c_hat, phi);
}
generated quantities {
vector[N] mu;
////vector[N] log_lik;
vector[N] y_rep;
mu = exp(c_hat);
for (i in 1:N) {
// log_lik[i] = neg_binomial_2_log_lpmf(y[i] | c_hat[i], phi);
y_rep[i] = neg_binomial_2_rng(c_hat[i], phi);
}
}
}