functions { /* * Alternative to poisson_log_rng() that avoids potential numerical problems during warmup */ int poisson_log_safe_rng(real eta) { real pois_rate = exp(eta); if (pois_rate >= exp(20.79)) return -9; return poisson_rng(pois_rate); } } data { // The input data are: int N; int complaints[N]; vector[N] traps; } parameters { // Our model accepts two parameters 'alpha' and 'beta' real alpha; real beta1; } model { // weakly informative priors: alpha ~ normal(log(4), 1); beta1 ~ normal(-0.25, 1); // likelihood: Consider poisson_log(eta), which is more efficient and stable alternative to poisson(exp(eta)) complaints ~ poisson_log(alpha + beta1 * traps); // Alternatively //complaints ~ poisson(exp(alpha + beta1 * traps)); } // Equivalently // model { // target += normal_lpdf(alpha|log(4),1); // target += normal_lpdf(beta1|-0.25,1); // target += poisson_log_lpmf(complaints|alpha + beta1 * traps); // } //Or // model { // beta1 ~ normal(-0.25, 1); // alpha ~ normal(log(4), 1); // target += poisson_log_lpmf(complaints|alpha + beta1 * traps); // } generated quantities { // sample predicted values from the model for posterior predictive checks int y_rep[N]; // for each obs we simulate a new value of y using the simulated parameters // the result will be a matrix with nrow = iter and ncol = n for (n in 1:N) { real eta_n = alpha + beta1 * traps[n]; y_rep[n] = poisson_log_safe_rng(eta_n); } }