Glad to hear a doc is coming for random.
In principle this should be possible by means of the "scaling property": Gamma(shape,scale) ~ scale*Gamma(shape,1)
FYI, you might recall that I started a mapping from math::distribution to random which I mention again because I've been refining it a bit. A Kolmogorov Smirnov accumulator is used for each distribution to verify that the mapping is correct for a particular parameter, for example, boost::mt19937 urng; typedef kolmovorov_smirnov::check_convergence<double> check_; math::gamma_distribution<double> d( 2, 3 ); check_()(6,10,10,d,make_random_generator(urng,d),os) prints gamma(2,3) (10,0.195231) (110,0.0720463) (1110,0.0190203) (11110,0.00706733) (111110,0.00400383) (1111110,0.000770538) https://svn.boost.org/svn/boost/sandbox/statistics/distribution_toolkit/ I've only dealt with what I needed but am open to suggestions.
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)
Some textbooks define Gamma(x|a,b) as prop to x^(a-1) exp(-b x) where b is called the "inverse scale" such as BDA by Andrew Gelman. Others, use x^(a-1) exp(- x/b) such as Wikipedia. I did not think through the implications here. It's just a mention in passing.