|
2329
|
1 REAL FUNCTION gennf(dfn,dfd,xnonc) |
|
|
2 |
|
|
3 C********************************************************************** |
|
|
4 C |
|
|
5 C REAL FUNCTION GENNF( DFN, DFD, XNONC ) |
|
|
6 C GENerate random deviate from the Noncentral F distribution |
|
|
7 C |
|
|
8 C |
|
|
9 C Function |
|
|
10 C |
|
|
11 C |
|
|
12 C Generates a random deviate from the noncentral F (variance ratio) |
|
|
13 C distribution with DFN degrees of freedom in the numerator, and DFD |
|
|
14 C degrees of freedom in the denominator, and noncentrality parameter |
|
|
15 C XNONC. |
|
|
16 C |
|
|
17 C |
|
|
18 C Arguments |
|
|
19 C |
|
|
20 C |
|
|
21 C DFN --> Numerator degrees of freedom |
|
|
22 C (Must be >= 1.0) |
|
|
23 C REAL DFN |
|
|
24 C DFD --> Denominator degrees of freedom |
|
|
25 C (Must be positive) |
|
|
26 C REAL DFD |
|
|
27 C |
|
|
28 C XNONC --> Noncentrality parameter |
|
|
29 C (Must be nonnegative) |
|
|
30 C REAL XNONC |
|
|
31 C |
|
|
32 C |
|
|
33 C Method |
|
|
34 C |
|
|
35 C |
|
|
36 C Directly generates ratio of noncentral numerator chisquare variate |
|
|
37 C to central denominator chisquare variate. |
|
|
38 C |
|
|
39 C********************************************************************** |
|
|
40 C .. Scalar Arguments .. |
|
|
41 REAL dfd,dfn,xnonc |
|
|
42 C .. |
|
|
43 C .. Local Scalars .. |
|
|
44 REAL xden,xnum |
|
|
45 LOGICAL qcond |
|
|
46 C .. |
|
|
47 C .. External Functions .. |
|
3188
|
48 C JJV changed the code to call SGAMMA and SNORM directly |
|
|
49 C REAL genchi,gennch |
|
|
50 C EXTERNAL genchi,gennch |
|
|
51 REAL sgamma,snorm |
|
|
52 EXTERNAL sgamma,snorm |
|
2329
|
53 C .. |
|
|
54 C .. Executable Statements .. |
|
3188
|
55 C JJV changed the argument checker to allow DFN = 1.0 |
|
|
56 C JJV in the same way as GENNCH was changed. |
|
|
57 qcond = dfn .LT. 1.0 .OR. dfd .LE. 0.0 .OR. xnonc .LT. 0.0 |
|
2329
|
58 IF (.NOT. (qcond)) GO TO 10 |
|
3188
|
59 WRITE (*,*) 'In GENNF - Either (1) Numerator DF < 1.0 or' |
|
|
60 WRITE (*,*) '(2) Denominator DF <= 0.0 or ' |
|
2329
|
61 WRITE (*,*) '(3) Noncentrality parameter < 0.0' |
|
|
62 WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd,'XNONC value: ', |
|
|
63 + xnonc |
|
3188
|
64 STOP 'Degrees of freedom or noncent param out of range in GENNF' |
|
|
65 |
|
|
66 C GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD ) |
|
|
67 C JJV changed this to call SGAMMA and SNORM directly |
|
|
68 C xnum = gennch(dfn,xnonc)/dfn |
|
|
69 10 IF (dfn.GE.1.000001) GO TO 20 |
|
|
70 C JJV case dfn = 1.0 - here I am treating dfn as exactly 1.0 |
|
|
71 xnum = (snorm() + sqrt(xnonc))**2 |
|
|
72 GO TO 30 |
|
2329
|
73 |
|
3188
|
74 C JJV case dfn > 1.0 |
|
|
75 20 xnum = (2.0*sgamma((dfn-1.0)/2.0) + (snorm()+sqrt(xnonc))**2)/dfn |
|
|
76 |
|
|
77 C xden = genchi(dfd)/dfd |
|
|
78 30 xden = 2.0*sgamma(dfd/2.0)/dfd |
|
|
79 |
|
|
80 C JJV changed constant so that it will not underflow at compile time |
|
|
81 C JJV while not slowing generator by using double precision or logs. |
|
|
82 C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 40 |
|
|
83 IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 40 |
|
2329
|
84 WRITE (*,*) ' GENNF - generated numbers would cause overflow' |
|
|
85 WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden |
|
3188
|
86 C JJV next 2 lines changed to maintain truncation of large deviates. |
|
|
87 C WRITE (*,*) ' GENNF returning 1.0E38' |
|
|
88 C gennf = 1.0E38 |
|
|
89 WRITE (*,*) ' GENNF returning 1.0E37' |
|
|
90 gennf = 1.0E37 |
|
|
91 GO TO 50 |
|
2329
|
92 |
|
3188
|
93 40 gennf = xnum/xden |
|
|
94 50 RETURN |
|
2329
|
95 |
|
|
96 END |
