GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
genf.f
Go to the documentation of this file.
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 ..
39 C JJV changed this code to call sgamma directly
40 C .. External Functions ..
41 C REAL genchi
42 C EXTERNAL genchi
43  REAL sgamma
44  EXTERNAL sgamma
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
50  CALL xstopx('Degrees of freedom nonpositive in GENF - abort!')
51
52  10 xnum = 2.0*sgamma(dfn/2.0)/dfn
53
54 C GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD )
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
60  WRITE (*,*) ' GENF - generated numbers would cause overflow'
61  WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden
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
67  go to 30
68
69  20 genf = xnum/xden
70  30 RETURN
71
72  END
real function genf(dfn, dfd)
Definition: genf.f:1