functions {
  int neg_binomial_2_log_safe_rng(real eta, real phi) {
    real gamma_rate = gamma_rng(phi, phi / exp(eta));
    if (gamma_rate >= exp(20.79))
      return -9;
    return poisson_rng(gamma_rate);
  }
}
data {
  int<lower=1> N;
  int<lower=0> complaints[N];
  vector<lower=0>[N] traps;
  vector<lower=0,upper=1>[N] super;
  vector<lower=0>[N] sqfoot;
}

transformed data{
    vector[N] log_sqfoot;
    log_sqfoot = log(sqfoot);
}

parameters {
  real alpha;
  real beta1;
  real beta2;
  real<lower=0> inv_phi;
}

transformed parameters {
  real phi = inv(inv_phi);
}

model {
  alpha ~ normal(log(4), 1);
  beta1 ~ normal(-0.25, 1);
  beta2 ~ normal(-0.5, 1);
  inv_phi ~ normal(0, 1);
  complaints ~ neg_binomial_2_log(alpha + beta1 * traps + beta2 * super + log_sqfoot, phi);
} 

generated quantities {
  int y_rep[N];
  vector[N] log_lik; 
  for (n in 1:N) {
    real eta_n = alpha + beta1 * traps[n] + beta2 * super[n] + log_sqfoot[n];
    y_rep[n] = neg_binomial_2_log_safe_rng(eta_n, phi);
    log_lik[n] = neg_binomial_2_log_lpmf(complaints[n]|eta_n, phi);
  }
}