SUBROUTINE GAMRAN(N,GAMMA,ISEED,X) C ******STILL NEEDS ALGORITHM WORK ****** C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE GAMMA DISTRIBUTION C WITH TAIL LENGTH PARAMETER VALUE = GAMMA. C THE PROTOTYPE GAMMA DISTRIBUTION USED C HEREIN HAS MEAN = GAMMA C AND STANDARD DEVIATION = SQRT(GAMMA). C THIS DISTRIBUTION IS DEFINED FOR ALL POSITIVE X, C AND HAS THE PROBABILITY DENSITY FUNCTION C F(X) = (1/CONSTANT) * (X**(GAMMA-1)) * EXP(-X) C WHERE THE CONSTANT = THE GAMMA FUNCTION EVALUATED C AT THE VALUE GAMMA. C INPUT ARGUMENTS--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. C --GAMMA = THE SINGLE PRECISION VALUE OF THE C TAIL LENGTH PARAMETER. C GAMMA SHOULD BE POSITIVE. C GAMMA SHOULD BE LARGER C THAN 1/3 (ALGORITHMIC RESTRICTION). C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE GAMMA DISTRIBUTION C WITH TAIL LENGTH PARAMETER VALUE = GAMMA. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C --GAMMA SHOULD BE POSITIVE. C --GAMMA SHOULD BE LARGER C THAN 1/3 (ALGORITHMIC RESTRICTION). C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN, NORRAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, EXP. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN (1977) C REFERENCES--GREENWOOD, 'A FAST GENERATOR FOR C GAMMA-DISTRIBUTED RANDOM VARIABLES', C COMPSTAT 1974, PROCEEDINGS IN C COMPUTATIONAL STATISTICS, VIENNA, C SEPTEMBER, 1974, PAGES 19-27. C --TOCHER, THE ART OF SIMULATION, C 1963, PAGES 24-27. C --HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS, C 1964, PAGES 36-37. C --WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY C PLOTS FOR THE GAMMA DISTRIBUTION', C TECHNOMETRICS, 1962, PAGES 1-15. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 166-206. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 68-73. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 952. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING DIVISION C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-3651 C NOTE--DATAPLOT IS A REGISTERED TRADEMARK C OF THE NATIONAL BUREAU OF STANDARDS. C THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED, C MODIFIED, OR OTHERWISE USED IN A CONTEXT C OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM. C LANGUAGE--ANSI FORTRAN (1966) C EXCEPTION--HOLLERITH STRINGS IN FORMAT STATEMENTS C DENOTED BY QUOTES RATHER THAN NH. C VERSION NUMBER--82/7 C ORIGINAL VERSION--NOVEMBER 1975. C UPDATED --FEBRUARY 1976. C UPDATED --JUNE 1978. C UPDATED --DECEMBER 1981. C UPDATED --MARCH 1982. C UPDATED --MAY 1982. C C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES------------------- C C--------------------------------------------------------------------- C DIMENSION X(*) C C--------------------------------------------------------------------- C CCCCC CHARACTER*4 IFEEDB CCCCC CHARACTER*4 IPRINT C CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW CCCCC COMMON /PRINT/IFEEDB,IPRINT C C-----DATA STATEMENTS------------------------------------------------- C DATA ATHIRD/0.3333333/ DATA SQRT3 /1.73205081/ C IPR=6 C C-----START POINT----------------------------------------------------- C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.1)GOTO50 IF(GAMMA.LE.0.0)GOTO60 IF(GAMMA.LE.0.33333333)GOTO65 GOTO90 50 WRITE(IPR, 5) WRITE(IPR,47)N RETURN 60 WRITE(IPR,15) WRITE(IPR,46)GAMMA RETURN 65 WRITE(IPR,16) WRITE(IPR,17) WRITE(IPR,46)GAMMA RETURN 90 CONTINUE 5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 GAMRAN SUBROUTINE IS NON-POSITIVE *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 GAMRAN SUBROUTINE IS NON-POSITIVE *****) 16 FORMAT(1H ,114H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 GAMRAN SUBROUTINE IS SMALLER THAN OR EQUAL TO 0.33333333 *****) 17 FORMAT(1H , 44H (ALGORITHMIC RESTIRCTION)) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C GENERATE N GAMMA DISTRIBUTION RANDOM NUMBERS C USING GREENWOOD'S REJECTION ALGORITHM-- C 1) GENERATE A NORMAL RANDOM NUMBER; C 2) TRANSFORM THE NORMAL VARIATE TO AN APPROXIMATE C GAMMA VARIATE USING THE WILSON-HILFERTY C APPROXIMATION (SEE THE JOHNSON AND KOTZ C REFERENCE, PAGE 176); C 3) FORM THE REJECTION FUNCTION VALUE, BASED C ON THE PROBABILITY DENSITY FUNCTION VALUE C OF THE ACTUAL DISTRIBUTION OF THE PSEUDO-GAMMA C VARIATE, AND THE PROBABILITY DENSITY FUNCTION VALUE C OF A TRUE GAMMA VARIATE. C 4) GENERATE A UNIFORM RANDOM NUMBER; C 5) IF THE UNIFORM RANDOM NUMBER IS LESS THAN C THE REJECTION FUNCTION VALUE, THEN ACCEPT C THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE; C IF THE UNIFORM RANDOM NUMBER IS LARGER THAN C THE REJECTION FUNCTION VALUE, THEN REJECT C THE PSEUDO-RANDOM NUMBER AS A GAMMA VARIATE. C A1=1.0/(9.0*GAMMA) B1=SQRT(A1) XN0=-SQRT3+B1 XG0=GAMMA*(1.0-A1+B1*XN0)**3 DO100I=1,N 150 CALL NORRAN(1,ISEED,XN) XG=GAMMA*(1.0-A1+B1*XN)**3 IF(XG.LT.0.0)GOTO150 TERM=(XG/XG0)**(GAMMA-ATHIRD) ARG=0.5*XN*XN-XG-0.5*XN0*XN0+XG0 FUNCT=TERM*EXP(ARG) CALL UNIRAN(1,ISEED,U) IF(U.LE.FUNCT)GOTO170 GOTO150 170 X(I)=XG 100 CONTINUE C RETURN END