Error evaluating the log probability at the initial value-Beverton Hold model works Ricker does not

Modeling
1

I am trying to fit stock-recruitment models using rstan: one beverton-holt, one ricker. The BH model runs “fine” (there are some divergences but I am just in inital fitting stages).

The Ricker model, is where I get the error: “Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: model42882f4835f2__namespace::log_prob: pRR[sym1__] is 6.78943e-05, but must be greater than or equal to 0.001 (in ‘string’, line 18, column 2 to column 27)”

I know this is because of how I restrict pRR in the model to be > 0.001, but I don’t understand why this is only happening in the Ricker model. How could the values of Alpha and Beta be so wrong that this error is occuring? I can’t share my data for privacy reasons. I am very new to STAN (~ 2 weeks of use).

My code for the BH is as follows:

BH<- "
data {
  int<lower=1> N; // number of datapoints
  int<lower=1> n_stream; // number of streams
  int<lower=1, upper=n_stream> stream[N]; //stream numeric index
  vector[N] SS; //spawners
  vector[N] lRR; //log recruits
}


parameters {
  real lAlpha[n_stream];
  real lBeta[n_stream];
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta;
  real<lower=0> sigma_obs;
  real lmu_alpha;
  real lmu_beta;

}
transformed parameters{
  real<lower=0.001> pRR[N];
  real lpRR[N];
  real<lower=0> Beta[n_stream];
  real<lower=0> Alpha[n_stream];

  Alpha = exp(lAlpha);
  Beta = exp(lBeta);

  for(i in 1:N){
    pRR[i] = exp(lAlpha[stream[i]])*SS[i]/(exp(lBeta[stream[i]])+SS[i]);
    lpRR[i] = log(pRR[i]);
  }
}
model {
  lmu_beta ~ normal(log(20),log(20));
  lmu_alpha ~ normal(log(500),log(500));
  //sigma_alpha ~ cauchy(0, 2.5);
  //sigma_y ~ cauchy(0, 2.5);
  sigma_obs ~ cauchy(0, 2.5);
  sigma_alpha ~ cauchy(0, 2.5);
  sigma_beta ~ cauchy(0, 2.5);

  lAlpha ~ normal(lmu_alpha, sigma_alpha);
  lBeta ~ normal(lmu_beta, sigma_beta);
  lRR ~ normal(lpRR,sigma_obs);
}
generated quantities{
  real<lower=0> pp_RR[N];
  for (i in 1:N) {
    pp_RR[i] = exp(normal_rng(lpRR[i] - 0.5 * sigma_obs^2, sigma_obs)); // generate posterior predictives
  }
}
"

dataIn <- list(N = nrow(dat),
               n_stream = max(dat$StreamNum),
               stream = dat$StreamNum,
               SS=dat$Redds.1km,
               lRR=log(dat$EST.1.1km))

#Fit model in parallel
fit_BH <- sampling(bhModel, 
                data=dataIn,
                thin=1,
                warmup=5000,
                iter=10000, 
                chains=4,
                control = list(max_treedepth = 15,
                               stepsize = 0.001,
                               adapt_delta = 0.99))

The stan code for Ricker is:

dataIn <- list(N = nrow(dat),
               n_stream = max(dat$StreamNum),
               stream = dat$StreamNum,
               SS=dat$Redds.1km,
               lRR=log(dat$EST.1.1km))

stan_Ricker_HG<- "

data {
  int<lower=1> N; // number of datapoints
  int<lower=1> n_stream; // number of streams
  int<lower=1, upper=n_stream> stream[N]; //stream numeric index
  vector[N] SS; //spawners
  vector[N] lRR; //log recruits
}

parameters {
  real lAlpha[n_stream];
  real lBeta[n_stream];
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta;
  real<lower=0> sigma_obs;
  real lmu_alpha;
  real lmu_beta;

}

transformed parameters{
  real<lower=0.001> pRR[N];
  real lpRR[N];
  real<upper=0> Beta[n_stream];
  real<lower=0> Alpha[n_stream];

  Alpha = exp(lAlpha);
  Beta = exp(lBeta);

  for(i in 1:N){
    pRR[i] = exp(Alpha[stream[i]])*SS[i]*exp(-1*Beta[stream[i]]*SS[i]);
    lpRR[i] = log(pRR[i]);
  }
}


model {
  lmu_beta ~ normal(log(20),log(20));
  lmu_alpha ~ normal(log(500),log(500));
  //sigma_alpha ~ cauchy(0, 2.5);
  //sigma_y ~ cauchy(0, 2.5);
  sigma_obs ~ cauchy(0, 2.5);
  sigma_alpha ~ cauchy(0, 2.5);
  sigma_beta ~ cauchy(0, 2.5);

  lAlpha ~ normal(lmu_alpha, sigma_alpha);
  lBeta ~ normal(lmu_beta, sigma_beta);
  lRR ~ normal(lpRR,sigma_obs);
}

generated quantities{
  real<lower=0> pp_RR[N];
  for (i in 1:N) {
    pp_RR[i] = exp(normal_rng(lpRR[i] - 0.5 * sigma_obs^2, sigma_obs)); // generate posterior predictives
  }
}
"

Ricker.fit_HG <- stan(model_code = stan_Ricker_HG, 
                      data=dataIn,
                      iter=10000, 
                      warmup=2000, 
                      thin=10, 
                      chains=4,
                      control=list(adapt_delta=0.9999,
                                   max_treedepth=20, 
                                   stepsize=0.1),cores=4)
2

The error occurs during initialization and Stan does not use the model until after initialzation. You should think of the first warmup draw as an arbitrary value that satisfies the constraints declared in the parameters block; in practice unconstrained real numbers initialize near zero (between -2 and 2). If your model needs an approximately correct starting point you need to set initial values explicitly.

By the way, the priors

  lmu_beta ~ normal(log(20),log(20));
  lmu_alpha ~ normal(log(500),log(500));

suggest you expect Beta ~ 20, Alpha ~ 500 but in the equation

    pRR[i] = exp(Alpha[stream[i]])*SS[i]*exp(-1*Beta[stream[i]]*SS[i]);

exp(500) is really, really large. In fact, too large for 64 bit floating point arithmetic, which is how Stan represents real numbers. Are you sure you didn’t mean to write

    pRR[i] = Alpha[stream[i]]*SS[i]*exp(-1*Beta[stream[i]]*SS[i]);

If you really do want exp(Alpha) you should write it as

    lpRR[i] = Alpha[stream[i]] + log(SS[i]) - Beta[stream[i]]*SS[i];
    pRR[i] = exp(lpRR[i]);

Regardless, probably the easiest way to adjust initial values is to declare offset; then the parameter initializes between offset-2 and offset+2 instead of around zero. For example

parameters {
  real<offset=log(500)> lAlpha[n_stream];
  real<offset=log(20)> lBeta[n_stream];
1 Like