Gamma probability distribution function.
Gamma(g,alpha,beta)
Generates a random variate following the Gamma probability distribution function, with parameters alpha and beta.
The pseudo-code algorithm used for this distribution is:
I. If alpha == 1, return y = exponential(1)
II. If alpha == 0.5, return y = normal(0,1)
III. If a < 1
        b = (e + alpha)/e
        1. Generate u1 as U(0,1). Let P = beta*u1. If P > 1, go to step 3. Otherwise, go to step 2.
        2. Let Y = P^(1/alpha), and generate u2 as U(0,1). If u2 &le e^(-Y), return y = Y. Otherwise, go to step 1.
        3. Let Y = -ln[(b - P)/alpha] and generate u2 as U(0,1). If u2 &le Y^(alpha-1), return y = Y. Otherwise, go to step 1.
IV. If a > 1
        a = 1/sqrt(2*alpha -1)
        b = alpha -ln(4)
        q = alpha + 1/a
        theta = 4.5
        d = 1 + ln(theta)
        1. Generate u1 and u2 as IID U(0,1).
        2. Let V = a*ln[u1/(1-u1)], Y = alpha*e^V, Z = u1^2*u2, and W = b + q*V - Y.
        3. If W + d - theta*Z &ge 0, return y = Y. Otherwise, go to step 4.
        3. If W &ge ln(Z), return y = Y. Otherwise, go to step 1.
V. x = beta*y;
(u,g) := Gamma(g,0.4,0.7);
function Gamma extends Var; end Gamma;