GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
gennf.f
Go to the documentation of this file.
1  REAL FUNCTION gennf(dfn,dfd,xnonc)
2 
3 C**********************************************************************
4 C
5 C REAL FUNCTION GENNF( DFN, DFD, XNONC )
6 C GENerate random deviate from the Noncentral F distribution
7 C
8 C
9 C Function
10 C
11 C
12 C Generates a random deviate from the noncentral F (variance ratio)
13 C distribution with DFN degrees of freedom in the numerator, and DFD
14 C degrees of freedom in the denominator, and noncentrality parameter
15 C XNONC.
16 C
17 C
18 C Arguments
19 C
20 C
21 C DFN --> Numerator degrees of freedom
22 C (Must be >= 1.0)
23 C REAL DFN
24 C DFD --> Denominator degrees of freedom
25 C (Must be positive)
26 C REAL DFD
27 C
28 C XNONC --> Noncentrality parameter
29 C (Must be nonnegative)
30 C REAL XNONC
31 C
32 C
33 C Method
34 C
35 C
36 C Directly generates ratio of noncentral numerator chisquare variate
37 C to central denominator chisquare variate.
38 C
39 C**********************************************************************
40 C .. Scalar Arguments ..
41  REAL dfd,dfn,xnonc
42 C ..
43 C .. Local Scalars ..
44  REAL xden,xnum
45  LOGICAL qcond
46 C ..
47 C .. External Functions ..
48 C JJV changed the code to call SGAMMA and SNORM directly
49 C REAL genchi,gennch
50 C EXTERNAL genchi,gennch
51  REAL sgamma,snorm
52  EXTERNAL sgamma,snorm
53 C ..
54 C .. Executable Statements ..
55 C JJV changed the argument checker to allow DFN = 1.0
56 C JJV in the same way as GENNCH was changed.
57  qcond = dfn .LT. 1.0 .OR. dfd .LE. 0.0 .OR. xnonc .LT. 0.0
58  IF (.NOT. (qcond)) go to 10
59  WRITE (*,*) 'In GENNF - Either (1) Numerator DF < 1.0 or'
60  WRITE (*,*) '(2) Denominator DF <= 0.0 or '
61  WRITE (*,*) '(3) Noncentrality parameter < 0.0'
62  WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd,'XNONC value: ',
63  + xnonc
64 
65  CALL xstopx
66  + ('Degrees of freedom or noncent param out of range in GENNF')
67 
68 C GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD )
69 C JJV changed this to call SGAMMA and SNORM directly
70 C xnum = gennch(dfn,xnonc)/dfn
71  10 IF (dfn.GE.1.000001) go to 20
72 C JJV case dfn = 1.0 - here I am treating dfn as exactly 1.0
73  xnum = (snorm() + sqrt(xnonc))**2
74  go to 30
75 
76 C JJV case dfn > 1.0
77  20 xnum = (2.0*sgamma((dfn-1.0)/2.0) + (snorm()+sqrt(xnonc))**2)/dfn
78 
79 C xden = genchi(dfd)/dfd
80  30 xden = 2.0*sgamma(dfd/2.0)/dfd
81 
82 C JJV changed constant so that it will not underflow at compile time
83 C JJV while not slowing generator by using double precision or logs.
84 C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 40
85  IF (.NOT. (xden.LE. (1.0e-37*xnum))) go to 40
86  WRITE (*,*) ' GENNF - generated numbers would cause overflow'
87  WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden
88 C JJV next 2 lines changed to maintain truncation of large deviates.
89 C WRITE (*,*) ' GENNF returning 1.0E38'
90 C gennf = 1.0E38
91  WRITE (*,*) ' GENNF returning 1.0E37'
92  gennf = 1.0e37
93  go to 50
94 
95  40 gennf = xnum/xden
96  50 RETURN
97 
98  END