00001 REAL FUNCTION genf(dfn,dfd) 00002 C********************************************************************** 00003 C 00004 C REAL FUNCTION GENF( DFN, DFD ) 00005 C GENerate random deviate from the F distribution 00006 C 00007 C 00008 C Function 00009 C 00010 C 00011 C Generates a random deviate from the F (variance ratio) 00012 C distribution with DFN degrees of freedom in the numerator 00013 C and DFD degrees of freedom in the denominator. 00014 C 00015 C 00016 C Arguments 00017 C 00018 C 00019 C DFN --> Numerator degrees of freedom 00020 C (Must be positive) 00021 C REAL DFN 00022 C DFD --> Denominator degrees of freedom 00023 C (Must be positive) 00024 C REAL DFD 00025 C 00026 C 00027 C Method 00028 C 00029 C 00030 C Directly generates ratio of chisquare variates 00031 C 00032 C********************************************************************** 00033 C .. Scalar Arguments .. 00034 REAL dfd,dfn 00035 C .. 00036 C .. Local Scalars .. 00037 REAL xden,xnum 00038 C .. 00039 C JJV changed this code to call sgamma directly 00040 C .. External Functions .. 00041 C REAL genchi 00042 C EXTERNAL genchi 00043 REAL sgamma 00044 EXTERNAL sgamma 00045 C .. 00046 C .. Executable Statements .. 00047 IF (.NOT. (dfn.LE.0.0.OR.dfd.LE.0.0)) GO TO 10 00048 WRITE (*,*) 'Degrees of freedom nonpositive in GENF - abort!' 00049 WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd 00050 CALL XSTOPX ('Degrees of freedom nonpositive in GENF - abort!') 00051 00052 10 xnum = 2.0*sgamma(dfn/2.0)/dfn 00053 00054 C GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD ) 00055 xden = 2.0*sgamma(dfd/2.0)/dfd 00056 C JJV changed constant so that it will not underflow at compile time 00057 C JJV while not slowing generator by using double precision or logs. 00058 C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 20 00059 IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 20 00060 WRITE (*,*) ' GENF - generated numbers would cause overflow' 00061 WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden 00062 C JJV next 2 lines changed to maintain truncation of large deviates. 00063 C WRITE (*,*) ' GENF returning 1.0E38' 00064 C genf = 1.0E38 00065 WRITE (*,*) ' GENF returning 1.0E37' 00066 genf = 1.0E37 00067 GO TO 30 00068 00069 20 genf = xnum/xden 00070 30 RETURN 00071 00072 END
1.7.1
