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