2329
|
1 REAL FUNCTION genf(dfn,dfd) |
|
2 C********************************************************************** |
|
3 C |
|
4 C REAL FUNCTION GENF( DFN, DFD ) |
|
5 C GENerate random deviate from the F distribution |
|
6 C |
|
7 C |
|
8 C Function |
|
9 C |
|
10 C |
|
11 C Generates a random deviate from the F (variance ratio) |
|
12 C distribution with DFN degrees of freedom in the numerator |
|
13 C and DFD degrees of freedom in the denominator. |
|
14 C |
|
15 C |
|
16 C Arguments |
|
17 C |
|
18 C |
|
19 C DFN --> Numerator degrees of freedom |
|
20 C (Must be positive) |
|
21 C REAL DFN |
|
22 C DFD --> Denominator degrees of freedom |
|
23 C (Must be positive) |
|
24 C REAL DFD |
|
25 C |
|
26 C |
|
27 C Method |
|
28 C |
|
29 C |
|
30 C Directly generates ratio of chisquare variates |
|
31 C |
|
32 C********************************************************************** |
|
33 C .. Scalar Arguments .. |
|
34 REAL dfd,dfn |
|
35 C .. |
|
36 C .. Local Scalars .. |
|
37 REAL xden,xnum |
|
38 C .. |
3188
|
39 C JJV changed this code to call sgamma directly |
2329
|
40 C .. External Functions .. |
3188
|
41 C REAL genchi |
|
42 C EXTERNAL genchi |
|
43 REAL sgamma |
|
44 EXTERNAL sgamma |
2329
|
45 C .. |
|
46 C .. Executable Statements .. |
|
47 IF (.NOT. (dfn.LE.0.0.OR.dfd.LE.0.0)) GO TO 10 |
|
48 WRITE (*,*) 'Degrees of freedom nonpositive in GENF - abort!' |
|
49 WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd |
3188
|
50 STOP 'Degrees of freedom nonpositive in GENF - abort!' |
|
51 |
|
52 10 xnum = 2.0*sgamma(dfn/2.0)/dfn |
2329
|
53 |
|
54 C GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD ) |
3188
|
55 xden = 2.0*sgamma(dfd/2.0)/dfd |
|
56 C JJV changed constant so that it will not underflow at compile time |
|
57 C JJV while not slowing generator by using double precision or logs. |
|
58 C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 20 |
|
59 IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 20 |
2329
|
60 WRITE (*,*) ' GENF - generated numbers would cause overflow' |
|
61 WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden |
3188
|
62 C JJV next 2 lines changed to maintain truncation of large deviates. |
|
63 C WRITE (*,*) ' GENF returning 1.0E38' |
|
64 C genf = 1.0E38 |
|
65 WRITE (*,*) ' GENF returning 1.0E37' |
|
66 genf = 1.0E37 |
2329
|
67 GO TO 30 |
|
68 |
|
69 20 genf = xnum/xden |
|
70 30 RETURN |
|
71 |
|
72 END |