00001 REAL FUNCTION gengam(a,r) 00002 C********************************************************************** 00003 C 00004 C REAL FUNCTION GENGAM( A, R ) 00005 C GENerates random deviates from GAMma distribution 00006 C 00007 C 00008 C Function 00009 C 00010 C 00011 C Generates random deviates from the gamma distribution whose 00012 C density is 00013 C (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X) 00014 C 00015 C 00016 C Arguments 00017 C 00018 C 00019 C JJV added the argument ranges supported 00020 C A --> Location parameter of Gamma distribution 00021 C REAL A ( A > 0 ) 00022 C 00023 C R --> Shape parameter of Gamma distribution 00024 C REAL R ( R > 0 ) 00025 C 00026 C 00027 C Method 00028 C 00029 C 00030 C Renames SGAMMA from TOMS as slightly modified by BWB to use RANF 00031 C instead of SUNIF. 00032 C 00033 C For details see: 00034 C (Case R >= 1.0) 00035 C Ahrens, J.H. and Dieter, U. 00036 C Generating Gamma Variates by a 00037 C Modified Rejection Technique. 00038 C Comm. ACM, 25,1 (Jan. 1982), 47 - 54. 00039 C Algorithm GD 00040 C 00041 C JJV altered the following to reflect sgamma argument ranges 00042 C (Case 0.0 < R < 1.0) 00043 C Ahrens, J.H. and Dieter, U. 00044 C Computer Methods for Sampling from Gamma, 00045 C Beta, Poisson and Binomial Distributions. 00046 C Computing, 12 (1974), 223-246/ 00047 C Adapted algorithm GS. 00048 C 00049 C********************************************************************** 00050 C .. Scalar Arguments .. 00051 REAL a,r 00052 C .. 00053 C .. External Functions .. 00054 REAL sgamma 00055 EXTERNAL sgamma 00056 C .. 00057 C .. Executable Statements .. 00058 00059 C JJV added argument value checker 00060 IF ( a.GT.0.0 .AND. r.GT.0.0 ) GO TO 10 00061 WRITE (*,*) 'In GENGAM - Either (1) Location param A <= 0.0 or' 00062 WRITE (*,*) '(2) Shape param R <= 0.0 - ABORT!' 00063 WRITE (*,*) 'A value: ',a,'R value: ',r 00064 CALL XSTOPX 00065 + ('Location or shape param out of range in GENGAM - ABORT!') 00066 C JJV end addition 00067 00068 10 gengam = sgamma(r)/a 00069 C gengam = gengam/a 00070 RETURN 00071 00072 END
1.7.1
