underdispersed Poisson alternatives

model-CMP.stan

functions {
  real CMP_log_lik(int Y, real mu, real nu) {
    return(nu * (Y * log(mu) - lgamma(Y+1)));  
  }
}

data {
  int C;
  int Y[C];
  int Count[C];
}

parameters {
  real<lower=0> mu;
  real<lower=0> nu;
}

model {
  real lp[11];
  for (c in 1:C)
    target += Count[c] * CMP_log_lik(Y[c], mu, nu);
  for (i in 0:10)
    lp[i+1] = CMP_log_lik(i, mu, nu);
  target += sum(Count) * (- log_sum_exp(lp));
}

model-DPO.stan

functions {
  real DPO_log_lik(int Y, real mu, real sigma) {
    real res;
    res = -0.5*log(sigma) - mu/sigma;
    if (Y == 0) {
      return(res);
    } else {
      return(res - Y + Y*log(Y) - lgamma(Y+1) + (Y/sigma)*(1 + log(mu) - log(Y)));  
    }
  }
}

data {
  int C;
  int Y[C];
  int Count[C];
}

parameters {
  real<lower=0> mu;
  real<lower=0> sigma;
}

model {
  real lp[11];
  for (c in 1:C)
    target += Count[c] * DPO_log_lik(Y[c], mu, sigma);
  for (i in 0:10)
    lp[i+1] = DPO_log_lik(i, mu, sigma);
  target += sum(Count) * (- log_sum_exp(lp));
}

model-LGP.stan

data {
  int C;
  int Y[C];
  int Count[C];
}

parameters {
  real<lower=0> theta;
  real<lower=-1, upper=1> lambda;
}

model {
  for (c in 1:C)
    target += Count[c] * (
      log(theta) + (Y[c]-1)*log(theta + lambda*Y[c]) - theta - lambda*Y[c] - lgamma(Y[c]+1)
    );
}

run-models.R

library(rstan)

d <- read.csv(file='z-data1.csv')
data <- list(C=nrow(d), Y=d$Y, Count=d$Count)

stanmodel1 <- stan_model(file='model-LGP.stan')
stanmodel2 <- stan_model(file='model-CMP.stan')
stanmodel3 <- stan_model(file='model-DPO.stan')

fit1 <- sampling(stanmodel1, data=data, seed=1234)
fit2 <- sampling(stanmodel2, data=data, seed=1234)
fit3 <- sampling(stanmodel3, data=data, seed=1234)

z-data1.csv