#lineson;
@A PROCEDURE FOR GENERATING GAMMA-DISTRIBUTED RANDOM NUMBERS
 USING VARIOUS METHODS DEPENDING ON THE VALUE OF THE ALPHA PARAMETER.
 WRITTEN BY RON MITTELHAMMER DECEMBER 1992 AND MODIFIED BY THOMAS
 HECKELEI JUNE 1993.

 PROVIDED WITHOUT GUARANTESS FOR PUBLIC NON-COMMERCIAL USE.

      N=NUMBER OF RANDOM NUMBERS DESIRED
      ALPHA=PARAMETER OF THE GAMMA FUNCTION
      BETA=PARAMETER IN EXPONENT OF E
@
PROC GAMGEN(ALPHA,BETA,N);
    LOCAL TEST, E,MULT,HOLD,NUMS,INVALPH,I,A,B,C,D,F,U,W,INVE,T,GX,
    MAXALPH,GAM;
    E=EXP(1);
    NUMS=ZEROS(N,1);
    TEST=TRUNC(ALPHA);
    IF ALPHA < 1; GOTO ONE;
    ELSEIF ALPHA == TEST AND ALPHA < 150; GOTO TWO;
    ELSE; GOTO THREE;
    ENDIF;

ONE:
     @USE THE INVERSE CDF METHOD HERE,TOGETHER WITH THE ACCEPTANCE/
      REJECTION METHOD. SEE AHRENS AND DIETER,COMPUTING,1974,P.223@

    INVALPH=1/ALPHA;
    INVE=1/E;
    MULT=E/(E+ALPHA);
    MAXALPH=1/(ALPHA*GAMMA(ALPHA)*MULT);
    I=1;
    DO UNTIL I > N;
        WIEDER:
        HOLD=RNDU(1,1);
        IF HOLD <= MULT;
            NUMS[I]=(HOLD/MULT)^INVALPH;
            GX=MULT*ALPHA*(NUMS[I]^(ALPHA-1));
        ELSE;
            NUMS[I]=-LN(INVE-((HOLD/MULT-1)*INVALPH));
            GX=MULT*ALPHA*EXP(-NUMS[I]);
        ENDIF;
        T=((1/(GAMMA(ALPHA)))*(NUMS[I]^(ALPHA-1))*EXP(-NUMS[I]))
          /(MAXALPH*GX);
        IF RNDU(1,1) > T; GOTO WIEDER; ENDIF;
        I=I+1;
    ENDO;
    RETP(NUMS*BETA);

TWO:

     @LOG OF PRODUCT OF ALPHA EXPONENTIAL DENSITIES FOR EVEN INTEGER
     ALPHA-VALUES. PROOF IS STRAIGHTFORWARD USING MOMENT GENERATING FUNCTION@

    GAM = -BETA*LN(PRODC(RNDU(ALPHA,N)));
    RETP(GAM);

THREE:
     @USING EITHER METHOD 1 OR METHOD 2 OF CHENG AND FEAST,APPLIED
      STATISTICS,1979,P.290, DEPENDING ON VALUE OF ALPHA.@

    A=ALPHA-1;
    B=(ALPHA-(1/(6*ALPHA)))/A;
    C=2/A;
    D=C+2;
    F=SQRT(ALPHA);
    I=1;
    DO UNTIL I > N;
        START:
        U=RNDU(2,1);
        IF ALPHA > 2.5;
            U[2]=U[1]+((1-1.86*U[2])/F);
        ENDIF;
        IF U[2] <= 0 OR U[2] >=1; GOTO START;
        ENDIF;
        W=B*U[1]/U[2];
        IF C*U[2]-D+W+(1/W) <=0; GOTO DOIT;
        ELSEIF C*LN(U[2])-LN(W)+W-1 >=0; GOTO START;
        ENDIF;
        DOIT:
        NUMS[I]=A*W;
        I=I+1;
    ENDO;
    RETP(NUMS*BETA);
ENDP;
END;

