rand_gamma

Generate a Gamma random variable
"Statistical Distributions", Evans, Hastings, Peacock, 2nd Edition,
Wiley, 1993, p.75-81

INPUTS:
(N,M) = size of array of random variables to generate
b = scale parameter > 0
c = shape parameter > 0

probability density function (pdf)
p(x) = (x/b)^(c-1) * exp(-x/b) / (b * gamma(c))

where gamma(c) is the gamma function (http://en.wikipedia.org/wiki/Gamma_function)

Basic stats of the gamma distribution
mean = b c
variance = b^2 c

generation method comes from
http://en.wikipedia.org/wiki/Gamma_distribution#Generating_gamma-distributed_random_variables
notation: theta = b
k = c
the algorithm exploits several properties of the gamma
(1) Gamma(b,1) ~ Exp(b) (an exponential random variable)
(2) sum_{i=1}^{n} Gamma(b,c_i) ~ Gamma(b,c) where c=sum_{i=1}^{n} c_i
plus acceptance-rejectance sampling