GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
gengam.f
Go to the documentation of this file.
1  REAL FUNCTION gengam(a,r)
2 C**********************************************************************
3 C
4 C REAL FUNCTION GENGAM( A, R )
5 C GENerates random deviates from GAMma distribution
6 C
7 C
8 C Function
9 C
10 C
11 C Generates random deviates from the gamma distribution whose
12 C density is
13 C (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X)
14 C
15 C
16 C Arguments
17 C
18 C
19 C JJV added the argument ranges supported
20 C A --> Location parameter of Gamma distribution
21 C REAL A ( A > 0 )
22 C
23 C R --> Shape parameter of Gamma distribution
24 C REAL R ( R > 0 )
25 C
26 C
27 C Method
28 C
29 C
30 C Renames SGAMMA from TOMS as slightly modified by BWB to use RANF
31 C instead of SUNIF.
32 C
33 C For details see:
34 C (Case R >= 1.0)
35 C Ahrens, J.H. and Dieter, U.
36 C Generating Gamma Variates by a
37 C Modified Rejection Technique.
38 C Comm. ACM, 25,1 (Jan. 1982), 47 - 54.
39 C Algorithm GD
40 C
41 C JJV altered the following to reflect sgamma argument ranges
42 C (Case 0.0 < R < 1.0)
43 C Ahrens, J.H. and Dieter, U.
44 C Computer Methods for Sampling from Gamma,
45 C Beta, Poisson and Binomial Distributions.
46 C Computing, 12 (1974), 223-246/
47 C Adapted algorithm GS.
48 C
49 C**********************************************************************
50 C .. Scalar Arguments ..
51  REAL a,r
52 C ..
53 C .. External Functions ..
54  REAL sgamma
55  EXTERNAL sgamma
56 C ..
57 C .. Executable Statements ..
58 
59 C JJV added argument value checker
60  IF ( a.GT.0.0 .AND. r.GT.0.0 ) go to 10
61  WRITE (*,*) 'In GENGAM - Either (1) Location param A <= 0.0 or'
62  WRITE (*,*) '(2) Shape param R <= 0.0 - ABORT!'
63  WRITE (*,*) 'A value: ',a,'R value: ',r
64  CALL xstopx
65  + ('Location or shape param out of range in GENGAM - ABORT!')
66 C JJV end addition
67 
68  10 gengam = sgamma(r)/a
69 C gengam = gengam/a
70  RETURN
71 
72  END