Gamma-count dist

count_gamma_reg.stan

functions {
  int sum0(vector x) {
    int len = num_elements(x);
    int res = 0;
    for (i in 1:len) if (x[i] == 0) res += 1;
    return res;
  }
}
data {
  int N;
  int p;
  matrix[N, p] X;
  vector<lower = 0>[N] y;
}
transformed data {
  int count_0 = sum0(y);
  int count_non0 = N - count_0;
  int idx_0[count_0];
  int idx_non0[count_non0];

  {
    int counter0 = 0;
    int counter_non0 = 0;
    for (i in 1:N) {
      if (y[i] == 0) {
        counter0 += 1;
        idx_0[counter0] = i;
      } else {
        counter_non0 += 1;
        idx_non0[counter_non0] = i;
      }
    }
  }
}
parameters {
  real ln_phi;
  real intercept;
  real<lower = 0> beta_sd;
  vector[p] beta;
}
model {
  real g_alpha = 1.0 / exp(ln_phi);
  vector[N] g_beta = 1.0 ./ exp(intercept + X * beta + ln_phi);

  ln_phi ~ student_t(3, 0, 5);
  intercept ~ student_t(3, 0, 5);
  beta_sd ~ student_t(3, 0, 5);
  beta ~ normal(0, beta_sd);

  // Winkelmann (1996) 10.1007/BF02926581
  target += sum(log1m(gamma_p(g_alpha, g_beta[idx_0])));
  target += sum(log(
    gamma_p(g_alpha * y[idx_non0], g_beta[idx_non0]) -
    gamma_p(g_alpha * y[idx_non0] + g_alpha, g_beta[idx_non0])));
  // // http://cursos.leg.ufpr.br/rmcd/models.html#gammacount
  // for (i in 1:N) {
  //   if (y[i] == 0) {
  //     target += log1m(gamma_cdf(1 | g_alpha, g_beta[i]));
  //   } else {
  //     target += log_diff_exp(
  //       gamma_lcdf(1 | y[i] * g_alpha, g_beta[i]),
  //       gamma_lcdf(1 | y[i] * g_alpha + g_alpha, g_beta[i]));
  //   }
  // }
}
generated quantities {
  real disp = exp(ln_phi);
}