* Bug #7296 fixed - Enabling %nan, %inf and -%inf for the cdf* functions
[scilab.git] / scilab / modules / statistics / src / dcdflib / cdfbin.f
1       SUBROUTINE cdfbin(which,p,q,s,xn,pr,ompr,status,bound)
2 C**********************************************************************
3 C
4 C      SUBROUTINE CDFBIN ( WHICH, P, Q, S, XN, PR, OMPR, STATUS, BOUND )
5 C               Cumulative Distribution Function
6 C                         BINomial distribution
7 C
8 C
9 C                              Function
10 C
11 C
12 C     Calculates any one parameter of the binomial
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 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 S,XN,PR and OMPR
23 C               iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
24 C               iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
25 C               iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
26 C                    INTEGER WHICH
27 C
28 C     P <--> The cumulation from 0 to S of the binomial distribution.
29 C            (Probablility of S or fewer successes in XN trials each
30 C            with probability of success PR.)
31 C            Input range: [0,1].
32 C                    DOUBLE PRECISION P
33 C
34 C     Q <--> 1-P.
35 C            Input range: [0, 1].
36 C            P + Q = 1.0.
37 C                    DOUBLE PRECISION Q
38 C
39 C     S <--> The number of successes observed.
40 C            Input range: [0, XN]
41 C            Search range: [0, XN]
42 C                    DOUBLE PRECISION S
43 C
44 C     XN  <--> The number of binomial trials.
45 C              Input range: (0, +infinity).
46 C              Search range: [1E-300, 1E300]
47 C                    DOUBLE PRECISION XN
48 C
49 C     PR  <--> The probability of success in each binomial trial.
50 C              Input range: [0,1].
51 C              Search range: [0,1]
52 C                    DOUBLE PRECISION PR
53 C
54 C     OMPR  <--> 1-PR
55 C              Input range: [0,1].
56 C              Search range: [0,1]
57 C              PR + OMPR = 1.0
58 C                    DOUBLE PRECISION OMPR
59 C
60 C     STATUS <-- 0 if calculation completed correctly
61 C               -I if input parameter number I is out of range
62 C                1 if answer appears to be lower than lowest
63 C                  search bound
64 C                2 if answer appears to be higher than greatest
65 C                  search bound
66 C                3 if P + Q .ne. 1
67 C                4 if PR + OMPR .ne. 1
68 C                    INTEGER STATUS
69 C
70 C     BOUND <-- Undefined if STATUS is 0
71 C
72 C               Bound exceeded by parameter number I if STATUS
73 C               is negative.
74 C
75 C               Lower search bound if STATUS is 1.
76 C
77 C               Upper search bound if STATUS is 2.
78 C
79 C
80 C                              Method
81 C
82 C
83 C     Formula  26.5.24    of   Abramowitz  and    Stegun,  Handbook   of
84 C     Mathematical   Functions (1966) is   used  to reduce the  binomial
85 C     distribution  to  the  cumulative incomplete    beta distribution.
86 C
87 C     Computation of other parameters involve a seach for a value that
88 C     produces  the desired  value  of P.   The search relies  on  the
89 C     monotinicity of P with the other parameter.
90 C
91 C
92 C**********************************************************************
93
94 C     .. Parameters ..
95       DOUBLE PRECISION atol
96       PARAMETER (atol=1.0D-50)
97       DOUBLE PRECISION tol
98       PARAMETER (tol=1.0D-13)
99       DOUBLE PRECISION zero,inf
100       PARAMETER (zero=1.0D-300,inf=1.0D300)
101       DOUBLE PRECISION one
102       PARAMETER (one=1.0D0)
103 C     ..
104 C     .. Scalar Arguments ..
105       DOUBLE PRECISION bound,p,q,pr,ompr,s,xn
106       INTEGER status,which
107 C     ..
108 C     .. Local Scalars ..
109       DOUBLE PRECISION fx,xhi,xlo,cum,ccum,pq,prompr
110       LOGICAL qhi,qleft,qporq
111 C     ..
112 C     .. External Functions ..
113       INTEGER vfinite
114       DOUBLE PRECISION spmpar
115       EXTERNAL spmpar
116 C     ..
117 C     .. External Subroutines ..
118       EXTERNAL dinvr,dstinv,dstzr,dzror,cumbin
119 C     ..
120 C     .. Executable Statements ..
121 C
122 C     Check arguments
123       IF ((which.GE.1).AND. (which.LE.4)) GO TO 30
124       IF (.NOT. (which.LT.1)) GO TO 10
125       bound = 1.0D0
126       GO TO 20
127
128    10 bound = 4.0D0
129    20 status = -1
130       RETURN
131
132    30 IF (which.EQ.1) GO TO 70
133 C
134 C     P
135 C
136       IF (ISANAN(p).EQ.1) THEN
137          CALL RETURNANANFORTRAN(pr)
138          CALL RETURNANANFORTRAN(ompr)
139          CALL RETURNANANFORTRAN(s)
140          CALL RETURNANANFORTRAN(xn)
141          RETURN
142       ENDIF
143       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
144       IF (.NOT. (p.LT.0.0D0)) GO TO 40
145       bound = 0.0D0
146       GO TO 50
147
148    40 bound = 1.0D0
149    50 status = -2
150       RETURN
151
152    60 CONTINUE
153    70 IF (which.EQ.1) GO TO 110
154 C
155 C     Q
156 C
157       IF (ISANAN(q).EQ.1) THEN
158          CALL RETURNANANFORTRAN(pr)
159          CALL RETURNANANFORTRAN(ompr)
160          CALL RETURNANANFORTRAN(s)
161          CALL RETURNANANFORTRAN(xn)
162          RETURN
163       ENDIF
164       IF (.NOT. ((q.LT.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
165       IF (.NOT. (q.LT.0.0D0)) GO TO 80
166       bound = 0.0D0
167       GO TO 90
168
169    80 bound = 1.0D0
170    90 status = -3
171       RETURN
172
173   100 CONTINUE
174   110 IF (which.EQ.3) GO TO 130
175 C
176 C     XN
177 C
178       IF (ISANAN(xn).EQ.1) THEN
179          CALL RETURNANANFORTRAN(p)
180          CALL RETURNANANFORTRAN(q)
181          CALL RETURNANANFORTRAN(s)
182          CALL RETURNANANFORTRAN(pr)
183          CALL RETURNANANFORTRAN(ompr)
184          RETURN
185       ENDIF
186       IF (vfinite(1,xn).EQ.0) xn = SIGN(inf,xn)
187       IF (.NOT. (xn.LE.0.0D0)) GO TO 120
188       bound = 0.0D0
189       status = -5
190       RETURN
191
192   120 CONTINUE
193   130 IF (which.EQ.2) GO TO 170
194 C
195 C     S
196 C
197       IF (ISANAN(s).EQ.1) THEN
198          CALL RETURNANANFORTRAN(p)
199          CALL RETURNANANFORTRAN(q)
200          CALL RETURNANANFORTRAN(xn)
201          CALL RETURNANANFORTRAN(pr)
202          CALL RETURNANANFORTRAN(ompr)
203          RETURN
204       ENDIF
205       IF (vfinite(1,s).EQ.0) s = SIGN(inf,s)
206       IF (.NOT. ((s.LT.0.0D0).OR. ((which.NE.3).AND.
207      +    (s.GT.xn)))) GO TO 160
208       IF (.NOT. (s.LT.0.0D0)) GO TO 140
209       bound = 0.0D0
210       GO TO 150
211
212   140 bound = xn
213   150 status = -4
214       RETURN
215
216   160 CONTINUE
217   170 IF (which.EQ.4) GO TO 210
218 C
219 C     PR
220 C
221       IF (ISANAN(pr).EQ.1) THEN
222          CALL RETURNANANFORTRAN(p)
223          CALL RETURNANANFORTRAN(q)
224          CALL RETURNANANFORTRAN(s)
225          CALL RETURNANANFORTRAN(xn)
226          RETURN
227       ENDIF
228       IF (.NOT. ((pr.LT.0.0D0).OR. (pr.GT.1.0D0))) GO TO 200
229       IF (.NOT. (pr.LT.0.0D0)) GO TO 180
230       bound = 0.0D0
231       GO TO 190
232
233   180 bound = 1.0D0
234   190 status = -6
235       RETURN
236
237   200 CONTINUE
238   210 IF (which.EQ.4) GO TO 250
239 C
240 C     OMPR
241 C
242       IF (ISANAN(ompr).EQ.1) THEN
243          CALL RETURNANANFORTRAN(p)
244          CALL RETURNANANFORTRAN(q)
245          CALL RETURNANANFORTRAN(s)
246          CALL RETURNANANFORTRAN(xn)
247          RETURN
248       ENDIF
249       IF (.NOT. ((ompr.LT.0.0D0).OR. (ompr.GT.1.0D0))) GO TO 240
250       IF (.NOT. (ompr.LT.0.0D0)) GO TO 220
251       bound = 0.0D0
252       GO TO 230
253
254   220 bound = 1.0D0
255   230 status = -7
256       RETURN
257
258   240 CONTINUE
259   250 IF (which.EQ.1) GO TO 290
260 C
261 C     P + Q
262 C
263       pq = p + q
264       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
265      +    (3.0D0*spmpar(1)))) GO TO 280
266       IF (.NOT. (pq.LT.0.0D0)) GO TO 260
267       bound = 0.0D0
268       GO TO 270
269
270   260 bound = 1.0D0
271   270 status = 3
272       RETURN
273
274   280 CONTINUE
275   290 IF (which.EQ.4) GO TO 330
276 C
277 C     PR + OMPR
278 C
279       prompr = pr + ompr
280       IF (.NOT. (abs(((prompr)-0.5D0)-0.5D0).GT.
281      +    (3.0D0*spmpar(1)))) GO TO 320
282       IF (.NOT. (prompr.LT.0.0D0)) GO TO 300
283       bound = 0.0D0
284       GO TO 310
285
286   300 bound = 1.0D0
287   310 status = 4
288       RETURN
289
290   320 CONTINUE
291   330 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
292 C
293 C     Select the minimum of P or Q
294 C
295 C
296 C     Calculate ANSWERS
297 C
298       IF ((1).EQ. (which)) THEN
299 C
300 C     Calculating P
301 C
302           CALL cumbin(s,xn,pr,ompr,p,q)
303           status = 0
304
305       ELSE IF ((2).EQ. (which)) THEN
306 C
307 C     Calculating S
308 C
309           s = 5.0D0
310           CALL dstinv(0.0D0,xn,0.5D0,0.5D0,5.0D0,atol,tol)
311           status = 0
312           CALL dinvr(status,s,fx,qleft,qhi)
313   340     IF (.NOT. (status.EQ.1)) GO TO 370
314           CALL cumbin(s,xn,pr,ompr,cum,ccum)
315           IF (.NOT. (qporq)) GO TO 350
316           fx = cum - p
317           GO TO 360
318
319   350     fx = ccum - q
320   360     CALL dinvr(status,s,fx,qleft,qhi)
321           GO TO 340
322
323   370     IF (.NOT. (status.EQ.-1)) GO TO 400
324           IF (.NOT. (qleft)) GO TO 380
325           status = 1
326           bound = 0.0D0
327           GO TO 390
328
329   380     status = 2
330           bound = xn
331   390     CONTINUE
332   400     CONTINUE
333
334       ELSE IF ((3).EQ. (which)) THEN
335 C
336 C     Calculating XN
337 C
338           xn = 5.0D0
339           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
340           status = 0
341           CALL dinvr(status,xn,fx,qleft,qhi)
342   410     IF (.NOT. (status.EQ.1)) GO TO 440
343           CALL cumbin(s,xn,pr,ompr,cum,ccum)
344           IF (.NOT. (qporq)) GO TO 420
345           fx = cum - p
346           GO TO 430
347
348   420     fx = ccum - q
349   430     CALL dinvr(status,xn,fx,qleft,qhi)
350           GO TO 410
351
352   440     IF (.NOT. (status.EQ.-1)) GO TO 470
353           IF (.NOT. (qleft)) GO TO 450
354           status = 1
355           bound = zero
356           GO TO 460
357
358   450     status = 2
359           bound = inf
360   460     CONTINUE
361   470     CONTINUE
362
363       ELSE IF ((4).EQ. (which)) THEN
364 C
365 C     Calculating PR and OMPR
366 C
367           CALL dstzr(0.0D0,1.0D0,atol,tol)
368           IF (.NOT. (qporq)) GO TO 500
369           status = 0
370           CALL dzror(status,pr,fx,xlo,xhi,qleft,qhi)
371           ompr = one - pr
372   480     IF (.NOT. (status.EQ.1)) GO TO 490
373           CALL cumbin(s,xn,pr,ompr,cum,ccum)
374           fx = cum - p
375           CALL dzror(status,pr,fx,xlo,xhi,qleft,qhi)
376           ompr = one - pr
377           GO TO 480
378
379   490     GO TO 530
380
381   500     status = 0
382           CALL dzror(status,ompr,fx,xlo,xhi,qleft,qhi)
383           pr = one - ompr
384   510     IF (.NOT. (status.EQ.1)) GO TO 520
385           CALL cumbin(s,xn,pr,ompr,cum,ccum)
386           fx = ccum - q
387           CALL dzror(status,ompr,fx,xlo,xhi,qleft,qhi)
388           pr = one - ompr
389           GO TO 510
390
391   520     CONTINUE
392   530     IF (.NOT. (status.EQ.-1)) GO TO 560
393           IF (.NOT. (qleft)) GO TO 540
394           status = 1
395           bound = 0.0D0
396           GO TO 550
397
398   540     status = 2
399           bound = 1.0D0
400   550     CONTINUE
401   560 END IF
402
403       RETURN
404
405       END