Generation of two-parameter Gamma variates

[Marco Guazzone <marco.guazzone@…>]

Actually, the code used to generate random Gamma variates cannot be used to generate variates for two-parameters Gamma distributions, that is for Gamma distributions with both shape and scale parameters (e.g., see ​http://en.wikipedia.org/wiki/Gamma_distribution or ​http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm).

In principle this should be possible by means of the scaling property: Gamma(shape,scale) ~ scale*Gamma(shape,1)

This would simply translate into

boost::mt19937 rng;
boost::gamma_distribution<> gamma(shape); // 1-parameter Gamma distribution
boost::variate_generator< boost::mt19937&, boost::gamma_distribution<> > rvg(rng, gamma);
double r = scale*rvg();

However, when shape == 1 the implementation uses the fact that Gamma(1) ~ Exp(1). This is rigth when scale==1 but not when scale != 1 since the exact relation is Gamma(1,scale) ~ Exp(1/scale)

[Marco Guazzone <marco.guazzone@…>]

R provides a possible implementation:

​https://svn.r-project.org/R/trunk/src/nmath/rgamma.c

[Steven Watanabe]

Might I point out that Exp(1/scale) ~ scale*Exp(1). In fact this is how Exp is implemented anyway. We probably lose slightly less accuracy by multiplying.

[Steven Watanabe]

Fixed in [63247].