GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
real function gennf(dfn, dfd, xnonc)
Definition: gennf.f:1
octave_value sqrt(void) const
Definition: ov.h:1200