00001 REAL FUNCTION gennf(dfn,dfd,xnonc) 00002 00003 C********************************************************************** 00004 C 00005 C REAL FUNCTION GENNF( DFN, DFD, XNONC ) 00006 C GENerate random deviate from the Noncentral F distribution 00007 C 00008 C 00009 C Function 00010 C 00011 C 00012 C Generates a random deviate from the noncentral F (variance ratio) 00013 C distribution with DFN degrees of freedom in the numerator, and DFD 00014 C degrees of freedom in the denominator, and noncentrality parameter 00015 C XNONC. 00016 C 00017 C 00018 C Arguments 00019 C 00020 C 00021 C DFN --> Numerator degrees of freedom 00022 C (Must be >= 1.0) 00023 C REAL DFN 00024 C DFD --> Denominator degrees of freedom 00025 C (Must be positive) 00026 C REAL DFD 00027 C 00028 C XNONC --> Noncentrality parameter 00029 C (Must be nonnegative) 00030 C REAL XNONC 00031 C 00032 C 00033 C Method 00034 C 00035 C 00036 C Directly generates ratio of noncentral numerator chisquare variate 00037 C to central denominator chisquare variate. 00038 C 00039 C********************************************************************** 00040 C .. Scalar Arguments .. 00041 REAL dfd,dfn,xnonc 00042 C .. 00043 C .. Local Scalars .. 00044 REAL xden,xnum 00045 LOGICAL qcond 00046 C .. 00047 C .. External Functions .. 00048 C JJV changed the code to call SGAMMA and SNORM directly 00049 C REAL genchi,gennch 00050 C EXTERNAL genchi,gennch 00051 REAL sgamma,snorm 00052 EXTERNAL sgamma,snorm 00053 C .. 00054 C .. Executable Statements .. 00055 C JJV changed the argument checker to allow DFN = 1.0 00056 C JJV in the same way as GENNCH was changed. 00057 qcond = dfn .LT. 1.0 .OR. dfd .LE. 0.0 .OR. xnonc .LT. 0.0 00058 IF (.NOT. (qcond)) GO TO 10 00059 WRITE (*,*) 'In GENNF - Either (1) Numerator DF < 1.0 or' 00060 WRITE (*,*) '(2) Denominator DF <= 0.0 or ' 00061 WRITE (*,*) '(3) Noncentrality parameter < 0.0' 00062 WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd,'XNONC value: ', 00063 + xnonc 00064 00065 CALL XSTOPX 00066 + ('Degrees of freedom or noncent param out of range in GENNF') 00067 00068 C GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD ) 00069 C JJV changed this to call SGAMMA and SNORM directly 00070 C xnum = gennch(dfn,xnonc)/dfn 00071 10 IF (dfn.GE.1.000001) GO TO 20 00072 C JJV case dfn = 1.0 - here I am treating dfn as exactly 1.0 00073 xnum = (snorm() + sqrt(xnonc))**2 00074 GO TO 30 00075 00076 C JJV case dfn > 1.0 00077 20 xnum = (2.0*sgamma((dfn-1.0)/2.0) + (snorm()+sqrt(xnonc))**2)/dfn 00078 00079 C xden = genchi(dfd)/dfd 00080 30 xden = 2.0*sgamma(dfd/2.0)/dfd 00081 00082 C JJV changed constant so that it will not underflow at compile time 00083 C JJV while not slowing generator by using double precision or logs. 00084 C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 40 00085 IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 40 00086 WRITE (*,*) ' GENNF - generated numbers would cause overflow' 00087 WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden 00088 C JJV next 2 lines changed to maintain truncation of large deviates. 00089 C WRITE (*,*) ' GENNF returning 1.0E38' 00090 C gennf = 1.0E38 00091 WRITE (*,*) ' GENNF returning 1.0E37' 00092 gennf = 1.0E37 00093 GO TO 50 00094 00095 40 gennf = xnum/xden 00096 50 RETURN 00097 00098 END
1.7.1
