Mercurial > hg > octave-lyh

annotate libcruft/ranlib/gennf.f @ 10535:3f973f6c841c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve sparse concatenation operator
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 20 Apr 2010 08:42:03 +0200
parents 0e71ead7359d
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
1 REAL FUNCTION gennf(dfn,dfd,xnonc)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
3 C**********************************************************************
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
4 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
5 C REAL FUNCTION GENNF( DFN, DFD, XNONC )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
6 C GENerate random deviate from the Noncentral F distribution
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
7 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
8 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
9 C Function
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
10 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
11 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
12 C Generates a random deviate from the noncentral F (variance ratio)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
13 C distribution with DFN degrees of freedom in the numerator, and DFD
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
14 C degrees of freedom in the denominator, and noncentrality parameter
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
15 C XNONC.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
16 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
17 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
18 C Arguments
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
19 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
20 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
21 C DFN --> Numerator degrees of freedom
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
22 C (Must be >= 1.0)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
23 C REAL DFN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
24 C DFD --> Denominator degrees of freedom
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
25 C (Must be positive)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
26 C REAL DFD
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
27 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
28 C XNONC --> Noncentrality parameter
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
29 C (Must be nonnegative)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
30 C REAL XNONC
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
31 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
32 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
33 C Method
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
34 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
35 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
36 C Directly generates ratio of noncentral numerator chisquare variate
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
37 C to central denominator chisquare variate.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
38 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
39 C**********************************************************************
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
40 C .. Scalar Arguments ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
41 REAL dfd,dfn,xnonc
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
42 C ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
43 C .. Local Scalars ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
44 REAL xden,xnum
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
45 LOGICAL qcond
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
46 C ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
47 C .. External Functions ..
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
48 C JJV changed the code to call SGAMMA and SNORM directly
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
49 C REAL genchi,gennch
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
50 C EXTERNAL genchi,gennch
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
51 REAL sgamma,snorm
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
52 EXTERNAL sgamma,snorm
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
53 C ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
54 C .. Executable Statements ..
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
55 C JJV changed the argument checker to allow DFN = 1.0
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
56 C JJV in the same way as GENNCH was changed.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
57 qcond = dfn .LT. 1.0 .OR. dfd .LE. 0.0 .OR. xnonc .LT. 0.0
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
58 IF (.NOT. (qcond)) GO TO 10
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
59 WRITE (*,*) 'In GENNF - Either (1) Numerator DF < 1.0 or'
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
60 WRITE (*,*) '(2) Denominator DF <= 0.0 or '
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
61 WRITE (*,*) '(3) Noncentrality parameter < 0.0'
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
62 WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd,'XNONC value: ',
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
63 + xnonc
10103
0e71ead7359d Use CALL XSTOPX instead of STOP in Fortran ranlib routines
Rik <rdrider0-list@yahoo.com>
parents: 3188
diff changeset
64
0e71ead7359d Use CALL XSTOPX instead of STOP in Fortran ranlib routines
Rik <rdrider0-list@yahoo.com>
parents: 3188
diff changeset
65 CALL XSTOPX
0e71ead7359d Use CALL XSTOPX instead of STOP in Fortran ranlib routines
Rik <rdrider0-list@yahoo.com>
parents: 3188
diff changeset
66 + ('Degrees of freedom or noncent param out of range in GENNF')
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
67
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
68 C GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD )
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
69 C JJV changed this to call SGAMMA and SNORM directly
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
70 C xnum = gennch(dfn,xnonc)/dfn
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
71 10 IF (dfn.GE.1.000001) GO TO 20
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
72 C JJV case dfn = 1.0 - here I am treating dfn as exactly 1.0
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
73 xnum = (snorm() + sqrt(xnonc))**2
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
74 GO TO 30
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
75
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
76 C JJV case dfn > 1.0
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
77 20 xnum = (2.0*sgamma((dfn-1.0)/2.0) + (snorm()+sqrt(xnonc))**2)/dfn
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
78
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
79 C xden = genchi(dfd)/dfd
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
80 30 xden = 2.0*sgamma(dfd/2.0)/dfd
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
81
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
82 C JJV changed constant so that it will not underflow at compile time
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
83 C JJV while not slowing generator by using double precision or logs.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
84 C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 40
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
85 IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 40
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
86 WRITE (*,*) ' GENNF - generated numbers would cause overflow'
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
87 WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
88 C JJV next 2 lines changed to maintain truncation of large deviates.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
89 C WRITE (*,*) ' GENNF returning 1.0E38'
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
90 C gennf = 1.0E38
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
91 WRITE (*,*) ' GENNF returning 1.0E37'
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
92 gennf = 1.0E37
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
93 GO TO 50
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
94
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
95 40 gennf = xnum/xden
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
96 50 RETURN
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
97
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
98 END