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