00001 REAL FUNCTION gennch(df,xnonc) 00002 C********************************************************************** 00003 C 00004 C REAL FUNCTION GENNCH( DF, XNONC ) 00005 C Generate random value of Noncentral CHIsquare variable 00006 C 00007 C 00008 C Function 00009 C 00010 C 00011 00012 C Generates random deviate from the distribution of a noncentral 00013 C chisquare with DF degrees of freedom and noncentrality parameter 00014 C XNONC. 00015 C 00016 C 00017 C Arguments 00018 C 00019 C 00020 C DF --> Degrees of freedom of the chisquare 00021 C (Must be >= 1.0) 00022 C REAL DF 00023 C 00024 C XNONC --> Noncentrality parameter of the chisquare 00025 C (Must be >= 0.0) 00026 C REAL XNONC 00027 C 00028 C 00029 C Method 00030 C 00031 C 00032 C Uses fact that noncentral chisquare is the sum of a chisquare 00033 C deviate with DF-1 degrees of freedom plus the square of a normal 00034 C deviate with mean sqrt(XNONC) and standard deviation 1. 00035 C 00036 C********************************************************************** 00037 C .. Scalar Arguments .. 00038 REAL df,xnonc 00039 C .. 00040 C .. External Functions .. 00041 C JJV changed these to call SGAMMA and SNORM directly 00042 C REAL genchi,gennor 00043 C EXTERNAL genchi,gennor 00044 REAL sgamma,snorm 00045 EXTERNAL sgamma,snorm 00046 C .. 00047 C .. Intrinsic Functions .. 00048 INTRINSIC sqrt 00049 C .. 00050 C JJV changed abort to df < 1, and added case: df = 1 00051 C .. Executable Statements .. 00052 IF (.NOT. (df.LT.1.0.OR.xnonc.LT.0.0)) GO TO 10 00053 WRITE (*,*) 'DF < 1 or XNONC < 0 in GENNCH - ABORT' 00054 WRITE (*,*) 'Value of DF: ',df,' Value of XNONC',xnonc 00055 CALL XSTOPX ('DF < 1 or XNONC < 0 in GENNCH - ABORT') 00056 00057 C JJV changed this to call SGAMMA and SNORM directly 00058 C gennch = genchi(df-1.0) + gennor(sqrt(xnonc),1.0)**2 00059 00060 10 IF (df.GE.1.000001) GO TO 20 00061 C JJV case DF = 1.0 00062 gennch = (snorm() + sqrt(xnonc))**2 00063 GO TO 30 00064 00065 C JJV case DF > 1.0 00066 20 gennch = 2.0*sgamma((df-1.0)/2.0) + (snorm() + sqrt(xnonc))**2 00067 30 RETURN 00068 00069 END
1.7.1
