* Bug #7296 fixed - Enabling %nan, %inf and -%inf for the cdf* functions
[scilab.git] / scilab / modules / statistics / src / dcdflib / cdff.f
1       SUBROUTINE cdff(which,p,q,f,dfn,dfd,status,bound)
2 C**********************************************************************
3 C
4 C      SUBROUTINE CDFF( WHICH, P, Q, F, DFN, DFD, STATUS, BOUND )
5 C               Cumulative Distribution Function
6 C               F distribution
7 C
8 C
9 C                              Function
10 C
11 C
12 C     Calculates any one parameter of the F distribution
13 C     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 four argument
20 C               values is to be calculated from the others.
21 C               Legal range: 1..4
22 C               iwhich = 1 : Calculate P and Q from F,DFN and DFD
23 C               iwhich = 2 : Calculate F from P,Q,DFN and DFD
24 C               iwhich = 3 : Calculate DFN from P,Q,F and DFD
25 C               iwhich = 4 : Calculate DFD from P,Q,F and DFN
26 C                    INTEGER WHICH
27 C
28 C       P <--> The integral from 0 to F of the f-density.
29 C              Input range: [0,1].
30 C                    DOUBLE PRECISION P
31 C
32 C       Q <--> 1-P.
33 C              Input range: (0, 1].
34 C              P + Q = 1.0.
35 C                    DOUBLE PRECISION Q
36 C
37 C       F <--> Upper limit of integration of the f-density.
38 C              Input range: [0, +infinity).
39 C              Search range: [0,1E300]
40 C                    DOUBLE PRECISION F
41 C
42 C     DFN < --> Degrees of freedom of the numerator sum of squares.
43 C               Input range: (0, +infinity).
44 C               Search range: [ 1E-300, 1E300]
45 C                    DOUBLE PRECISION DFN
46 C
47 C     DFD < --> Degrees of freedom of the denominator sum of squares.
48 C               Input range: (0, +infinity).
49 C               Search range: [ 1E-300, 1E300]
50 C                    DOUBLE PRECISION DFD
51 C
52 C     STATUS <-- 0 if calculation completed correctly
53 C               -I if input parameter number I is out of range
54 C                1 if answer appears to be lower than lowest
55 C                  search bound
56 C                2 if answer appears to be higher than greatest
57 C                  search bound
58 C                3 if P + Q .ne. 1
59 C                    INTEGER STATUS
60 C
61 C     BOUND <-- Undefined if STATUS is 0
62 C
63 C               Bound exceeded by parameter number I if STATUS
64 C               is negative.
65 C
66 C               Lower search bound if STATUS is 1.
67 C
68 C               Upper search bound if STATUS is 2.
69 C
70 C
71 C                              Method
72 C
73 C
74 C     Formula   26.6.2   of   Abramowitz   and   Stegun,  Handbook  of
75 C     Mathematical  Functions (1966) is used to reduce the computation
76 C     of the  cumulative  distribution function for the  F  variate to
77 C     that of an incomplete beta.
78 C
79 C     Computation of other parameters involve a seach for a value that
80 C     produces  the desired  value  of P.   The search relies  on  the
81 C     monotinicity of P with the other parameter.
82 C
83 C                              WARNING
84 C
85 C     The value of the  cumulative  F distribution is  not necessarily
86 C     monotone in  either degrees of freedom.  There  thus may  be two
87 C     values  that  provide a given CDF  value.   This routine assumes
88 C     monotonicity and will find an arbitrary one of the two values.
89 C
90 C**********************************************************************
91 C     .. Parameters ..
92       DOUBLE PRECISION tol
93       PARAMETER (tol=1.0D-13)
94       DOUBLE PRECISION atol
95       PARAMETER (atol=1.0D-50)
96       DOUBLE PRECISION zero,inf
97       PARAMETER (zero=1.0D-300,inf=1.0D300)
98 C     ..
99 C     .. Scalar Arguments ..
100       DOUBLE PRECISION bound,dfd,dfn,f,p,q
101       INTEGER status,which
102 C     ..
103 C     .. Local Scalars ..
104       DOUBLE PRECISION pq,fx,cum,ccum
105       LOGICAL qhi,qleft,qporq
106 C     ..
107 C     .. External Functions ..
108       INTEGER vfinite
109       DOUBLE PRECISION spmpar
110       EXTERNAL spmpar
111 C     ..
112 C     .. External Subroutines ..
113       EXTERNAL dinvr,dstinv,cumf
114 C     ..
115 C     .. Executable Statements ..
116 C
117 C     Check arguments
118 C
119       IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
120       IF (.NOT. (which.LT.1)) GO TO 10
121       bound = 1.0D0
122       GO TO 20
123
124    10 bound = 4.0D0
125    20 status = -1
126       RETURN
127
128    30 IF (which.EQ.1) GO TO 70
129 C
130 C     P
131 C
132       IF (ISANAN(p).EQ.1) THEN
133          CALL RETURNANANFORTRAN(f)
134          CALL RETURNANANFORTRAN(dfn)
135          CALL RETURNANANFORTRAN(dfd)
136          RETURN
137       ENDIF
138       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
139       IF (.NOT. (p.LT.0.0D0)) GO TO 40
140       bound = 0.0D0
141       GO TO 50
142
143    40 bound = 1.0D0
144    50 status = -2
145       RETURN
146
147    60 CONTINUE
148    70 IF (which.EQ.1) GO TO 110
149 C
150 C     Q
151 C
152       IF (ISANAN(q).EQ.1) THEN
153          CALL RETURNANANFORTRAN(f)
154          CALL RETURNANANFORTRAN(dfn)
155          CALL RETURNANANFORTRAN(dfd)
156          RETURN
157       ENDIF
158       IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
159       IF (.NOT. (q.LE.0.0D0)) GO TO 80
160       bound = 0.0D0
161       GO TO 90
162
163    80 bound = 1.0D0
164    90 status = -3
165       RETURN
166
167   100 CONTINUE
168   110 IF (which.EQ.2) GO TO 130
169 C
170 C     F
171 C
172       IF (ISANAN(f).EQ.1) THEN
173          CALL RETURNANANFORTRAN(p)
174          CALL RETURNANANFORTRAN(q)
175          CALL RETURNANANFORTRAN(dfn)
176          CALL RETURNANANFORTRAN(dfd)
177          RETURN
178       ENDIF
179       IF (vfinite(1,f).EQ.0) then
180          IF (which.EQ.1) then
181             IF (f.GT.0) then
182                p = 1
183                q = 0
184                RETURN
185             ENDIF
186          ELSE
187             f = SIGN(1D300,f)
188          ENDIF
189       ENDIF
190       IF (.NOT. (f.LT.0.0D0)) GO TO 120
191       bound = 0.0D0
192       status = -4
193       RETURN
194
195   120 CONTINUE
196   130 IF (which.EQ.3) GO TO 150
197 C
198 C     DFN
199 C
200       IF (ISANAN(dfn).EQ.1) THEN
201          CALL RETURNANANFORTRAN(p)
202          CALL RETURNANANFORTRAN(q)
203          CALL RETURNANANFORTRAN(f)
204          CALL RETURNANANFORTRAN(dfd)
205          RETURN
206       ENDIF
207       IF (vfinite(1,dfn).EQ.0) dfn = SIGN(inf,dfn)
208       IF (.NOT. (dfn.LE.0.0D0)) GO TO 140
209       bound = 0.0D0
210       status = -5
211       RETURN
212
213   140 CONTINUE
214   150 IF (which.EQ.4) GO TO 170
215 C
216 C     DFD
217 C
218       IF (ISANAN(dfd).EQ.1) THEN
219          CALL RETURNANANFORTRAN(p)
220          CALL RETURNANANFORTRAN(q)
221          CALL RETURNANANFORTRAN(f)
222          CALL RETURNANANFORTRAN(dfn)
223          RETURN
224       ENDIF
225       IF (vfinite(1,dfd).EQ.0) dfd = SIGN(inf,dfd)
226       IF (.NOT. (dfd.LE.0.0D0)) GO TO 160
227       bound = 0.0D0
228       status = -6
229       RETURN
230
231   160 CONTINUE
232   170 IF (which.EQ.1) GO TO 210
233 C
234 C     P + Q
235 C
236       pq = p + q
237       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
238      +    (3.0D0*spmpar(1)))) GO TO 200
239       IF (.NOT. (pq.LT.0.0D0)) GO TO 180
240       bound = 0.0D0
241       GO TO 190
242
243   180 bound = 1.0D0
244   190 status = 3
245       RETURN
246
247   200 CONTINUE
248   210 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
249 C
250 C     Select the minimum of P or Q
251 C
252 C
253 C     Calculate ANSWERS
254 C
255       IF ((1).EQ. (which)) THEN
256 C
257 C     Calculating P
258 C
259           CALL cumf(f,dfn,dfd,p,q)
260           status = 0
261
262       ELSE IF ((2).EQ. (which)) THEN
263 C
264 C     Calculating F
265 C
266           f = 5.0D0
267           CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
268           status = 0
269           CALL dinvr(status,f,fx,qleft,qhi)
270   220     IF (.NOT. (status.EQ.1)) GO TO 250
271           CALL cumf(f,dfn,dfd,cum,ccum)
272           IF (.NOT. (qporq)) GO TO 230
273           fx = cum - p
274           GO TO 240
275
276   230     fx = ccum - q
277   240     CALL dinvr(status,f,fx,qleft,qhi)
278           GO TO 220
279
280   250     IF (.NOT. (status.EQ.-1)) GO TO 280
281           IF (.NOT. (qleft)) GO TO 260
282           status = 1
283           bound = 0.0D0
284           GO TO 270
285
286   260     status = 2
287           bound = inf
288   270     CONTINUE
289   280     CONTINUE
290
291       ELSE IF ((3).EQ. (which)) THEN
292 C
293 C     Calculating DFN
294 C
295           dfn = 5.0D0
296           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
297           status = 0
298           CALL dinvr(status,dfn,fx,qleft,qhi)
299   290     IF (.NOT. (status.EQ.1)) GO TO 320
300           CALL cumf(f,dfn,dfd,cum,ccum)
301           IF (.NOT. (qporq)) GO TO 300
302           fx = cum - p
303           GO TO 310
304
305   300     fx = ccum - q
306   310     CALL dinvr(status,dfn,fx,qleft,qhi)
307           GO TO 290
308
309   320     IF (.NOT. (status.EQ.-1)) GO TO 350
310           IF (.NOT. (qleft)) GO TO 330
311           status = 1
312           bound = zero
313           GO TO 340
314
315   330     status = 2
316           bound = inf
317   340     CONTINUE
318   350     CONTINUE
319
320       ELSE IF ((4).EQ. (which)) THEN
321 C
322 C     Calculating DFD
323 C
324           dfd = 5.0D0
325           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
326           status = 0
327           CALL dinvr(status,dfd,fx,qleft,qhi)
328   360     IF (.NOT. (status.EQ.1)) GO TO 390
329           CALL cumf(f,dfn,dfd,cum,ccum)
330           IF (.NOT. (qporq)) GO TO 370
331           fx = cum - p
332           GO TO 380
333
334   370     fx = ccum - q
335   380     CALL dinvr(status,dfd,fx,qleft,qhi)
336           GO TO 360
337
338   390     IF (.NOT. (status.EQ.-1)) GO TO 420
339           IF (.NOT. (qleft)) GO TO 400
340           status = 1
341           bound = zero
342           GO TO 410
343
344   400     status = 2
345           bound = inf
346   410     CONTINUE
347   420 END IF
348
349       RETURN
350
351       END