Bug? Poisson rng is generating negative numbers?

Modeling
#1

This question is a follow-up to my previous question. The following stan code is a minimum reproducible example of the problem that I am encountering:

data{
  int<lower=0> N;
  int<lower=0> offset[N];
  int<lower=0> y[N];
  int<lower=0> pred_N;
  int<lower=0> pred_offset[pred_N];
}

parameters {
  real alpha_poisson ;                
}


model{
  alpha_poisson ~ normal(0,1) ; 
  for(n in 1:N) target += poisson_log_lpmf(y[n] | alpha_poisson + log(offset[n]));
}

generated quantities {
  int<lower=0>  y_sim[N] ;
  int<lower=0>  y_pred[N] ;
  
  for(n in 1:N){
    if(alpha_poisson + log(offset[n])<20.7944)
	      y_sim[n]=poisson_log_rng(alpha_poisson + log(offset[n]));
	     else
	      y_sim[n]=999999;
  }
  
  for(n in 1:pred_N){
	     if(alpha_poisson + log(pred_offset[n])<20.7944)
	      y_pred[n]=poisson_log_rng(alpha_poisson + log(pred_offset[n]));
	     else
	      y_pred[n]=999999; 
  }
}

When I run this I get the following error:

Chain 1: Exception: model2130306f11b_poisson_namespace::write_array: y_pred[k0__] is -2147483648, but must be greater than or equal to 0  (in 'model2130306f11b_poisson' at line 21)

Am I doing something wrong or is this a bug in the random number generator?

Thanks a lot for the help!

Posterior predictive checking
#2

It overflowed.

1 Like
#3

Thanks for responding so quickly @bgoodri . Is this a problem with my code then? do you have any suggestion for how to solve this problem?

#4

I would first try it with save_warmup = FALSE and hope that this problem only occurs during the warmup phase, in which case you don’t really need to do anything about it.

#5

I tried running fit <- sampling(stan_model, data = stan_data, seed=9782, save_warmup = FALSE). Alas, I’m still getting the same error :_(

#6

OK, well the real issue is that your Poisson intensities are getting much too big. So, I would figure that out first.

#7

perhaps setting the comparison to something definitely smaller than 20.7944 - (30 log(2)) - how about 20.7 ?

1 Like
#8

Setting it to 20.7 did allows me to run the code! Now i need to understand why am I getting numbers sooo big that make absolute no sense in my context.

#9

your problem is on line 21 -

int<lower=0> y_pred[N] ;

size of y_pred should be N_pred

all values of y_pred have been initialized to std::numeric_limits<int>::min(). it seems that N_pred is less than N, so some values of y_pred fail to meet the constraint.
when N_pred is less than or equal to N you don’t have this problem.

1 Like
#10

Thanks a lot @mitzimorris and @bgoodri! In case this is helpful for someone else in the future, this code is working:

data{
  int<lower=0> N;
  int<lower=0> offset[N];
  int<lower=0> y[N];
  int<lower=0> pred_N;
  int<lower=0> pred_offset[pred_N];
}

parameters {
  real alpha_poisson ;                
}


model{
  alpha_poisson ~ normal(0,1) ; 
  for(n in 1:N) target += poisson_log_lpmf(y[n] | alpha_poisson + log(offset[n]));
}

generated quantities {
  int<lower=0>  y_sim[N] ;
  int<lower=0>  y_pred[pred_N] ;
  
  for(n in 1:N){
    if(alpha_poisson + log(offset[n])<20.7)
	      y_sim[n]=poisson_log_rng(alpha_poisson + log(offset[n]));
	     else
	      y_sim[n]=999999;
  }
  
  for(n in 1:pred_N){
	     if(alpha_poisson + log(pred_offset[n])<20.7)
	      y_pred[n]=poisson_log_rng(alpha_poisson + log(pred_offset[n]));
	     else
	      y_pred[n]=999999;
  }
}