1       SUBROUTINE cdffnc(which,p,q,f,dfn,dfd,phonc,status,bound)
2 C**********************************************************************
3 C
4 C      SUBROUTINE CDFFNC( WHICH, P, Q, F, DFN, DFD, PNONC, STATUS, BOUND
5 C               Cumulative Distribution Function
6 C               Non-central F distribution
7 C
8 C
9 C                              Function
10 C
11 C
12 C     Calculates any one parameter of the Non-central F
13 C     distribution given values for the others.
14 C
15 C
16 C                              Arguments
17 C
18 C
19 C     WHICH --> Integer indicating which of the next five argument
20 C               values is to be calculated from the others.
21 C               Legal range: 1..5
22 C               iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC
23 C               iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC
24 C               iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC
25 C               iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC
26 C               iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD
27 C                    INTEGER WHICH
28 C
29 C       P <--> The integral from 0 to F of the non-central f-density.
30 C              Input range: [0,1-1E-16).
31 C                    DOUBLE PRECISION P
32 C
33 C       Q <--> 1-P.
34 C            Q is not used by this subroutine and is only included
35 C            for similarity with other cdf* routines.
36 C                    DOUBLE PRECISION Q
37 C
38 C       F <--> Upper limit of integration of the non-central f-density.
39 C              Input range: [0, +infinity).
40 C              Search range: [0,1E300]
41 C                    DOUBLE PRECISION F
42 C
43 C     DFN < --> Degrees of freedom of the numerator sum of squares.
44 C               Input range: (0, +infinity).
45 C               Search range: [ 1E-300, 1E300]
46 C                    DOUBLE PRECISION DFN
47 C
48 C     DFD < --> Degrees of freedom of the denominator sum of squares.
49 C               Must be in range: (0, +infinity).
50 C               Input range: (0, +infinity).
51 C               Search range: [ 1E-300, 1E300]
52 C                    DOUBLE PRECISION DFD
53 C
54 C     PNONC <-> The non-centrality parameter
55 C               Input range: [0,infinity)
56 C               Search range: [0,1E4]
57 C                    DOUBLE PRECISION PHONC
58 C
59 C     STATUS <-- 0 if calculation completed correctly
60 C               -I if input parameter number I is out of range
61 C                1 if answer appears to be lower than lowest
62 C                  search bound
63 C                2 if answer appears to be higher than greatest
64 C                  search bound
65 C                3 if P + Q .ne. 1
66 C                    INTEGER STATUS
67 C
68 C     BOUND <-- Undefined if STATUS is 0
69 C
70 C               Bound exceeded by parameter number I if STATUS
71 C               is negative.
72 C
73 C               Lower search bound if STATUS is 1.
74 C
75 C               Upper search bound if STATUS is 2.
76 C
77 C
78 C                              Method
79 C
80 C
81 C     Formula  26.6.20   of   Abramowitz   and   Stegun,  Handbook  of
82 C     Mathematical  Functions (1966) is used to compute the cumulative
83 C     distribution function.
84 C
85 C     Computation of other parameters involve a seach for a value that
86 C     produces  the desired  value  of P.   The search relies  on  the
87 C     monotinicity of P with the other parameter.
88 C
89 C                            WARNING
90 C
91 C     The computation time  required for this  routine is proportional
92 C     to the noncentrality  parameter  (PNONC).  Very large  values of
93 C     this parameter can consume immense  computer resources.  This is
94 C     why the search range is bounded by 10,000.
95 C
96 C                              WARNING
97 C
98 C     The  value  of the  cumulative  noncentral F distribution is not
99 C     necessarily monotone in either degrees  of freedom.  There  thus
100 C     may be two values that provide a given  CDF value.  This routine
101 C     assumes monotonicity  and will find  an arbitrary one of the two
102 C     values.
103 C
104 C**********************************************************************
105 C     .. Parameters ..
106       DOUBLE PRECISION tent4
107       PARAMETER (tent4=1.0D4)
108       DOUBLE PRECISION tol
109       PARAMETER (tol=1.0D-13)
110       DOUBLE PRECISION atol
111       PARAMETER (atol=1.0D-50)
112       DOUBLE PRECISION zero,one,inf
113       PARAMETER (zero=1.0D-300,one=1.0D0-1.0D-16,inf=1.0D300)
114 C     ..
115 C     .. Scalar Arguments ..
116       DOUBLE PRECISION bound,dfd,dfn,f,p,q,phonc
117       INTEGER status,which
118 C     ..
119 C     .. Local Scalars ..
120       DOUBLE PRECISION fx,cum,ccum
121       LOGICAL qhi,qleft
122 C     ..
123 C     .. External Subroutines ..
124       EXTERNAL dinvr,dstinv,cumfnc
125 C     ..
126 C     .. Executable Statements ..
127 C
128 C     Check arguments
129 C
130       IF (.NOT. ((which.LT.1).OR. (which.GT.5))) GO TO 30
131       IF (.NOT. (which.LT.1)) GO TO 10
132       bound = 1.0D0
133       GO TO 20
135    10 bound = 5.0D0
136    20 status = -1
137       RETURN
139    30 IF (which.EQ.1) GO TO 70
140 C
141 C     P
142 C
143       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.one))) GO TO 60
144       IF (.NOT. (p.LT.0.0D0)) GO TO 40
145       bound = 0.0D0
146       GO TO 50
148    40 bound = one
149    50 status = -2
150       RETURN
152    60 CONTINUE
153    70 IF (which.EQ.2) GO TO 90
154 C
155 C     F
156 C
157       IF (.NOT. (f.LT.0.0D0)) GO TO 80
158       bound = 0.0D0
159       status = -4
160       RETURN
162    80 CONTINUE
163    90 IF (which.EQ.3) GO TO 110
164 C
165 C     DFN
166 C
167       IF (.NOT. (dfn.LE.0.0D0)) GO TO 100
168       bound = 0.0D0
169       status = -5
170       RETURN
172   100 CONTINUE
173   110 IF (which.EQ.4) GO TO 130
174 C
175 C     DFD
176 C
177       IF (.NOT. (dfd.LE.0.0D0)) GO TO 120
178       bound = 0.0D0
179       status = -6
180       RETURN
182   120 CONTINUE
183   130 IF (which.EQ.5) GO TO 150
184 C
185 C     PHONC
186 C
187       IF (.NOT. (phonc.LT.0.0D0)) GO TO 140
188       bound = 0.0D0
189       status = -7
190       RETURN
192   140 CONTINUE
193 C
195 C
196   150 IF ((1).EQ. (which)) THEN
197 C
198 C     Calculating P
199 C
200           CALL cumfnc(f,dfn,dfd,phonc,p,q)
201           status = 0
203       ELSE IF ((2).EQ. (which)) THEN
204 C
205 C     Calculating F
206 C
207           f = 5.0D0
208           CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
209           status = 0
210           CALL dinvr(status,f,fx,qleft,qhi)
211   160     IF (.NOT. (status.EQ.1)) GO TO 170
212           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
213           fx = cum - p
214           CALL dinvr(status,f,fx,qleft,qhi)
215           GO TO 160
217   170     IF (.NOT. (status.EQ.-1)) GO TO 200
218           IF (.NOT. (qleft)) GO TO 180
219           status = 1
220           bound = 0.0D0
221           GO TO 190
223   180     status = 2
224           bound = inf
225   190     CONTINUE
226   200     CONTINUE
228       ELSE IF ((3).EQ. (which)) THEN
229 C
230 C     Calculating DFN
231 C
232           dfn = 5.0D0
233           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
234           status = 0
235           CALL dinvr(status,dfn,fx,qleft,qhi)
236   210     IF (.NOT. (status.EQ.1)) GO TO 220
237           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
238           fx = cum - p
239           CALL dinvr(status,dfn,fx,qleft,qhi)
240           GO TO 210
242   220     IF (.NOT. (status.EQ.-1)) GO TO 250
243           IF (.NOT. (qleft)) GO TO 230
244           status = 1
245           bound = zero
246           GO TO 240
248   230     status = 2
249           bound = inf
250   240     CONTINUE
251   250     CONTINUE
253       ELSE IF ((4).EQ. (which)) THEN
254 C
255 C     Calculating DFD
256 C
257           dfd = 5.0D0
258           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
259           status = 0
260           CALL dinvr(status,dfd,fx,qleft,qhi)
261   260     IF (.NOT. (status.EQ.1)) GO TO 270
262           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
263           fx = cum - p
264           CALL dinvr(status,dfd,fx,qleft,qhi)
265           GO TO 260
267   270     IF (.NOT. (status.EQ.-1)) GO TO 300
268           IF (.NOT. (qleft)) GO TO 280
269           status = 1
270           bound = zero
271           GO TO 290
273   280     status = 2
274           bound = inf
275   290     CONTINUE
276   300     CONTINUE
278       ELSE IF ((5).EQ. (which)) THEN
279 C
280 C     Calculating PHONC
281 C
282           phonc = 5.0D0
283           CALL dstinv(0.0D0,tent4,0.5D0,0.5D0,5.0D0,atol,tol)
284           status = 0
285           CALL dinvr(status,phonc,fx,qleft,qhi)
286   310     IF (.NOT. (status.EQ.1)) GO TO 320
287           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
288           fx = cum - p
289           CALL dinvr(status,phonc,fx,qleft,qhi)
290           GO TO 310
292   320     IF (.NOT. (status.EQ.-1)) GO TO 350
293           IF (.NOT. (qleft)) GO TO 330
294           status = 1
295           bound = 0.0D0
296           GO TO 340
298   330     status = 2
299           bound = tent4
300   340     CONTINUE
301   350 END IF
303       RETURN
305       END