/*
Gamma distributed random number generator


Procedure: Gamma distributed random numbers
Written by Ronald Schoenberg for public non-commercial use, 14 July 1995.

Below is a Gauss proc for generating gamma distributed random variables.
It is a transcription of a Fortran program appearing in _Principles_of_
_Random_Variate_Generation_ by John Dagpunar.
*/

************************************************************************
*    Ronald     *      Aptech Systems        *   rons@Aptech.com       *
*  Schoenberg   * University of Washington   *   rons@u.washington.edu *
************************************************************************


/*
**  rndgam.src      Gamma distributed random numbers
**
**  Transcribed from Fortran that appeared in _Principles_of_Random_
**  _Variate_Generation_ by John Dagpunar, Clarendon Press, 1988.
**
**  The source code below is in the public domain, and the author
**  cannot be held responsible for correctness or any actual use.
**
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**
**  FORMAT
**          x = RNDGAM(R,C,alpha);
**
**  INPUT
**          R - number of rows of matrix of random numbers
**
**          C - number of columns of matrix of random numbers
**
**          alpha - scalar, Rx1 column vector, 1xC row vector, or RxC
**                  matrix, scale parameters
**
**  OUTPUT
**          x - RxC matrix of gamma distributed random numbers
**
**
**  REMARKS
**
**   If alpha is
**       1x1 scalar, used for all variates
**       Rx1 vector, row of alpha used for each row of result
**       1xC vector, column of alpha used for each column of result
**       RxC matrix, element of alpha used for each element of result
**
**
**   MOMENTS:
**
**   E(x) = VAR(x) = alpha
**
**
**   TWO PARAMETER MODEL:
**
**   For two parameter gamma random variables, i.e., where
**
**           theta^alpha * x ^ (alpha-1) * exp(-x / theta)
**    p(x) = ----------------------------------------------
**                theta^alpha * gamma(alpha)
**
**   generate a random variate using RNDG() with alpha and
**   multiply the result by theta.  Thus
**
**      E(x) = alpha * theta
**    VAR(x) = alpha * theta^2
**
**   For example,
**
**       theta = 1.2;
**       alpha = .8;
**       x = theta.*rndg(4,4,alpha);
**
**   creates a 4x4 matrix of gamma distributed variables with
**   theta = 1.2 and alpha = .8.
**
*/


proc rndgam(r,c,alpha);
    local i,j,z,alpha0;

    if rows(alpha) == c and cols(alpha) == 1;
        alpha = alpha';
    endif;
    if cols(alpha) == r and rows(alpha) == 1;
        alpha = alpha';
    endif;

    z = zeros(r,c);
    i = 1;
    do until i > r;
        j = 1;
        do until j > c;
            if cols(alpha) == c and rows(alpha) == r;
                alpha0 = alpha[i,j];
            elseif rows(alpha) == r;
                alpha0 = alpha[i,.];
            elseif cols(alpha) == c;
                alpha0 = alpha[.,j];
            else;
                alpha0 = alpha[1,1];
            endif;
            z[i,j] = _rndgm(alpha0);
            j = j + 1;
       endo;
       i = i + 1;
    endo;
    retp(z);
endp;



proc _rndgm(alpha);
    if alpha > 1;
        retp(_rndgm0(alpha));
    elseif alpha > 0 and alpha < 1;
        retp(_rndgm1(alpha));
    elseif alpha == 1;
        retp(_rndgm2);
    else;
        retp(error(0));
    endif;
endp;


proc _rndgm0(a);
    local b,h,r,g,d,f,x;
    b = a - 1;
    h = sqrt(3*a-.75);

A1:

    r = rndu(1,1);
    g = r - r * r;
    if g <= 0;
       goto A1;
    endif;
    f = (r - .5) * h / sqrt(g);
    x = b + f;
    if x <= 0;
        goto A1;
    endif;
    d = 64 * g * (rndu(1,1) * g) ^ 2;
    if d <= 0 or d * x < x - 2*f*f;
        retp(x);
    endif;
    if ln(d) >= 2*(b*ln(x/b)-f);
        goto A1;
    endif;
    retp(x);
endp;


proc _rndgm1(a);
   local b,p,c,uf,u,s,d,r,x,w;

   u = __macheps;
   b = 1 - a;
   p = b / (b + a * exp(-b));
   if p == 1;
     retp(error(0));
   endif;
   c = 1 / a;
   uf = p * (u / b) ^ a;
   d = b * ln(b);

B1:

   r = rndu(1,1);
   if r + u >= 1;
       goto B1;
   elseif r > p;
       x = b - ln((1 - r) / (1 - p));
       w = b * ln(x) - d;
   elseif r > uf;
       x = b * (r / p) ^ c;
       w = x;
   else;
       retp(error(0));
   endif;

   r = rndu(1,1);
   if (1 - r <= w) and (r > 0) and (r*(w+1) >= 1 or -ln(r) <= w);
       goto B1;
   endif;
   retp(x);
endp;


proc _rndgm2;
   local r;
   r = rndu(1,1);
   do while r <= 0;
       r = rndu(1,1);
   endo;
   retp(-ln(r));
endp;
