2dd86ac9b51e4add680fad4835d030f35458879c
[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       DOUBLE PRECISION spmpar
123       EXTERNAL spmpar
124 C     ..
125 C     .. External Subroutines ..
126       EXTERNAL dinvr,dstinv,dstzr,dzror,cumnbn
127 C     ..
128 C     .. Executable Statements ..
129 C
130 C     Check arguments
131 C
132       IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
133       IF (.NOT. (which.LT.1)) GO TO 10
134       bound = 1.0D0
135       GO TO 20
136
137    10 bound = 4.0D0
138    20 status = -1
139       RETURN
140
141    30 IF (which.EQ.1) GO TO 70
142 C
143 C     P
144 C
145       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
146       IF (.NOT. (p.LT.0.0D0)) GO TO 40
147       bound = 0.0D0
148       GO TO 50
149
150    40 bound = 1.0D0
151    50 status = -2
152       RETURN
153
154    60 CONTINUE
155    70 IF (which.EQ.1) GO TO 110
156 C
157 C     Q
158 C
159       IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
160       IF (.NOT. (q.LE.0.0D0)) GO TO 80
161       bound = 0.0D0
162       GO TO 90
163
164    80 bound = 1.0D0
165    90 status = -3
166       RETURN
167
168   100 CONTINUE
169   110 IF (which.EQ.2) GO TO 130
170 C
171 C     S
172 C
173       IF (.NOT. (s.LT.0.0D0)) GO TO 120
174       bound = 0.0D0
175       status = -4
176       RETURN
177
178   120 CONTINUE
179   130 IF (which.EQ.3) GO TO 150
180 C
181 C     XN
182 C
183       IF (.NOT. (xn.LT.0.0D0)) GO TO 140
184       bound = 0.0D0
185       status = -5
186       RETURN
187
188   140 CONTINUE
189   150 IF (which.EQ.4) GO TO 190
190 C
191 C     PR
192 C
193       IF (.NOT. ((pr.LT.0.0D0).OR. (pr.GT.1.0D0))) GO TO 180
194       IF (.NOT. (pr.LT.0.0D0)) GO TO 160
195       bound = 0.0D0
196       GO TO 170
197
198   160 bound = 1.0D0
199   170 status = -6
200       RETURN
201
202   180 CONTINUE
203   190 IF (which.EQ.4) GO TO 230
204 C
205 C     OMPR
206 C
207       IF (.NOT. ((ompr.LT.0.0D0).OR. (ompr.GT.1.0D0))) GO TO 220
208       IF (.NOT. (ompr.LT.0.0D0)) GO TO 200
209       bound = 0.0D0
210       GO TO 210
211
212   200 bound = 1.0D0
213   210 status = -7
214       RETURN
215
216   220 CONTINUE
217   230 IF (which.EQ.1) GO TO 270
218 C
219 C     P + Q
220 C
221       pq = p + q
222       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
223      +    (3.0D0*spmpar(1)))) GO TO 260
224       IF (.NOT. (pq.LT.0.0D0)) GO TO 240
225       bound = 0.0D0
226       GO TO 250
227
228   240 bound = 1.0D0
229   250 status = 3
230       RETURN
231
232   260 CONTINUE
233   270 IF (which.EQ.4) GO TO 310
234 C
235 C     PR + OMPR
236 C
237       prompr = pr + ompr
238       IF (.NOT. (abs(((prompr)-0.5D0)-0.5D0).GT.
239      +    (3.0D0*spmpar(1)))) GO TO 300
240       IF (.NOT. (prompr.LT.0.0D0)) GO TO 280
241       bound = 0.0D0
242       GO TO 290
243
244   280 bound = 1.0D0
245   290 status = 4
246       RETURN
247
248   300 CONTINUE
249   310 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
250 C
251 C     Select the minimum of P or Q
252 C
253 C
254 C     Calculate ANSWERS
255 C
256       IF ((1).EQ. (which)) THEN
257 C
258 C     Calculating P
259 C
260           CALL cumnbn(s,xn,pr,ompr,p,q)
261           status = 0
262
263       ELSE IF ((2).EQ. (which)) THEN
264 C
265 C     Calculating S
266 C
267           s = 5.0D0
268           CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
269           status = 0
270           CALL dinvr(status,s,fx,qleft,qhi)
271   320     IF (.NOT. (status.EQ.1)) GO TO 350
272           CALL cumnbn(s,xn,pr,ompr,cum,ccum)
273           IF (.NOT. (qporq)) GO TO 330
274           fx = cum - p
275           GO TO 340
276
277   330     fx = ccum - q
278   340     CALL dinvr(status,s,fx,qleft,qhi)
279           GO TO 320
280
281   350     IF (.NOT. (status.EQ.-1)) GO TO 380
282           IF (.NOT. (qleft)) GO TO 360
283           status = 1
284           bound = 0.0D0
285           GO TO 370
286
287   360     status = 2
288           bound = inf
289   370     CONTINUE
290   380     CONTINUE
291
292       ELSE IF ((3).EQ. (which)) THEN
293 C
294 C     Calculating XN
295 C
296           xn = 5.0D0
297           CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
298           status = 0
299           CALL dinvr(status,xn,fx,qleft,qhi)
300   390     IF (.NOT. (status.EQ.1)) GO TO 420
301           CALL cumnbn(s,xn,pr,ompr,cum,ccum)
302           IF (.NOT. (qporq)) GO TO 400
303           fx = cum - p
304           GO TO 410
305
306   400     fx = ccum - q
307   410     CALL dinvr(status,xn,fx,qleft,qhi)
308           GO TO 390
309
310   420     IF (.NOT. (status.EQ.-1)) GO TO 450
311           IF (.NOT. (qleft)) GO TO 430
312           status = 1
313           bound = 0.0D0
314           GO TO 440
315
316   430     status = 2
317           bound = inf
318   440     CONTINUE
319   450     CONTINUE
320
321       ELSE IF ((4).EQ. (which)) THEN
322 C
323 C     Calculating PR and OMPR
324 C
325           CALL dstzr(0.0D0,1.0D0,atol,tol)
326           IF (.NOT. (qporq)) GO TO 480
327           status = 0
328           CALL dzror(status,pr,fx,xlo,xhi,qleft,qhi)
329           ompr = one - pr
330   460     IF (.NOT. (status.EQ.1)) GO TO 470
331           CALL cumnbn(s,xn,pr,ompr,cum,ccum)
332           fx = cum - p
333           CALL dzror(status,pr,fx,xlo,xhi,qleft,qhi)
334           ompr = one - pr
335           GO TO 460
336
337   470     GO TO 510
338
339   480     status = 0
340           CALL dzror(status,ompr,fx,xlo,xhi,qleft,qhi)
341           pr = one - ompr
342   490     IF (.NOT. (status.EQ.1)) GO TO 500
343           CALL cumnbn(s,xn,pr,ompr,cum,ccum)
344           fx = ccum - q
345           CALL dzror(status,ompr,fx,xlo,xhi,qleft,qhi)
346           pr = one - ompr
347           GO TO 490
348
349   500     CONTINUE
350   510     IF (.NOT. (status.EQ.-1)) GO TO 540
351           IF (.NOT. (qleft)) GO TO 520
352           status = 1
353           bound = 0.0D0
354           GO TO 530
355
356   520     status = 2
357           bound = 1.0D0
358   530     CONTINUE
359   540 END IF
360
361       RETURN
362
363       END