90435cd76ea443c2da8ca998d4e8ced2bd9a41e7
[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       DOUBLE PRECISION spmpar
109       EXTERNAL spmpar
110 C     ..
111 C     .. External Subroutines ..
112       EXTERNAL dinvr,dstinv,cumf
113 C     ..
114 C     .. Executable Statements ..
115 C
116 C     Check arguments
117 C
118       IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
119       IF (.NOT. (which.LT.1)) GO TO 10
120       bound = 1.0D0
121       GO TO 20
122
123    10 bound = 4.0D0
124    20 status = -1
125       RETURN
126
127    30 IF (which.EQ.1) GO TO 70
128 C
129 C     P
130 C
131       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
132       IF (.NOT. (p.LT.0.0D0)) GO TO 40
133       bound = 0.0D0
134       GO TO 50
135
136    40 bound = 1.0D0
137    50 status = -2
138       RETURN
139
140    60 CONTINUE
141    70 IF (which.EQ.1) GO TO 110
142 C
143 C     Q
144 C
145       IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
146       IF (.NOT. (q.LE.0.0D0)) GO TO 80
147       bound = 0.0D0
148       GO TO 90
149
150    80 bound = 1.0D0
151    90 status = -3
152       RETURN
153
154   100 CONTINUE
155   110 IF (which.EQ.2) GO TO 130
156 C
157 C     F
158 C
159       IF (.NOT. (f.LT.0.0D0)) GO TO 120
160       bound = 0.0D0
161       status = -4
162       RETURN
163
164   120 CONTINUE
165   130 IF (which.EQ.3) GO TO 150
166 C
167 C     DFN
168 C
169       IF (.NOT. (dfn.LE.0.0D0)) GO TO 140
170       bound = 0.0D0
171       status = -5
172       RETURN
173
174   140 CONTINUE
175   150 IF (which.EQ.4) GO TO 170
176 C
177 C     DFD
178 C
179       IF (.NOT. (dfd.LE.0.0D0)) GO TO 160
180       bound = 0.0D0
181       status = -6
182       RETURN
183
184   160 CONTINUE
185   170 IF (which.EQ.1) GO TO 210
186 C
187 C     P + Q
188 C
189       pq = p + q
190       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
191      +    (3.0D0*spmpar(1)))) GO TO 200
192       IF (.NOT. (pq.LT.0.0D0)) GO TO 180
193       bound = 0.0D0
194       GO TO 190
195
196   180 bound = 1.0D0
197   190 status = 3
198       RETURN
199
200   200 CONTINUE
201   210 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
202 C
203 C     Select the minimum of P or Q
204 C
205 C
206 C     Calculate ANSWERS
207 C
208       IF ((1).EQ. (which)) THEN
209 C
210 C     Calculating P
211 C
212           CALL cumf(f,dfn,dfd,p,q)
213           status = 0
214
215       ELSE IF ((2).EQ. (which)) THEN
216 C
217 C     Calculating F
218 C
219           f = 5.0D0
220           CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
221           status = 0
222           CALL dinvr(status,f,fx,qleft,qhi)
223   220     IF (.NOT. (status.EQ.1)) GO TO 250
224           CALL cumf(f,dfn,dfd,cum,ccum)
225           IF (.NOT. (qporq)) GO TO 230
226           fx = cum - p
227           GO TO 240
228
229   230     fx = ccum - q
230   240     CALL dinvr(status,f,fx,qleft,qhi)
231           GO TO 220
232
233   250     IF (.NOT. (status.EQ.-1)) GO TO 280
234           IF (.NOT. (qleft)) GO TO 260
235           status = 1
236           bound = 0.0D0
237           GO TO 270
238
239   260     status = 2
240           bound = inf
241   270     CONTINUE
242   280     CONTINUE
243
244       ELSE IF ((3).EQ. (which)) THEN
245 C
246 C     Calculating DFN
247 C
248           dfn = 5.0D0
249           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
250           status = 0
251           CALL dinvr(status,dfn,fx,qleft,qhi)
252   290     IF (.NOT. (status.EQ.1)) GO TO 320
253           CALL cumf(f,dfn,dfd,cum,ccum)
254           IF (.NOT. (qporq)) GO TO 300
255           fx = cum - p
256           GO TO 310
257
258   300     fx = ccum - q
259   310     CALL dinvr(status,dfn,fx,qleft,qhi)
260           GO TO 290
261
262   320     IF (.NOT. (status.EQ.-1)) GO TO 350
263           IF (.NOT. (qleft)) GO TO 330
264           status = 1
265           bound = zero
266           GO TO 340
267
268   330     status = 2
269           bound = inf
270   340     CONTINUE
271   350     CONTINUE
272
273       ELSE IF ((4).EQ. (which)) THEN
274 C
275 C     Calculating DFD
276 C
277           dfd = 5.0D0
278           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
279           status = 0
280           CALL dinvr(status,dfd,fx,qleft,qhi)
281   360     IF (.NOT. (status.EQ.1)) GO TO 390
282           CALL cumf(f,dfn,dfd,cum,ccum)
283           IF (.NOT. (qporq)) GO TO 370
284           fx = cum - p
285           GO TO 380
286
287   370     fx = ccum - q
288   380     CALL dinvr(status,dfd,fx,qleft,qhi)
289           GO TO 360
290
291   390     IF (.NOT. (status.EQ.-1)) GO TO 420
292           IF (.NOT. (qleft)) GO TO 400
293           status = 1
294           bound = zero
295           GO TO 410
296
297   400     status = 2
298           bound = inf
299   410     CONTINUE
300   420 END IF
301
302       RETURN
303
304       END