Negative binomial simulation

[rachaelmountain]

In the below, I think it should be beta=p/(1-p). The rand_gamma() function uses arma::randg to simulate from the gamma distribution. beta=(1-p)/p would be correct if we were defining beta as the rate parameter of the gamma, but randg defines beta as the scale parameter (=1/rate) - this is the opposite of the default for rgamma in R which takes beta as the rate. Could someone check this?

int rand_NegBin(double rate, double dispersion)
{
  if (dispersion != 1)
  {
    double size=1/dispersion;
    double p=size/(size+rate);
    double alpha=size;
    double beta=(1-p)/p;

    arma::vec lambda_arma = rand_gamma(alpha,beta);
    double lambda = lambda_arma(0);

    int x=rand_Poisson(lambda);
    return(x);
  }

  if (dispersion == 1)
    {
    int x=rand_Poisson(rate);
    return(x);
    }

  return(0);
}