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**********************************************************************
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
128    10 bound = 4.0D0
129    20 status = -1
130       RETURN
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
148    40 bound = 1.0D0
149    50 status = -2
150       RETURN
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
169    80 bound = 1.0D0
170    90 status = -3
171       RETURN
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
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
212   140 bound = xn
213   150 status = -4
214       RETURN
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
233   180 bound = 1.0D0
234   190 status = -6
235       RETURN
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
254   220 bound = 1.0D0
255   230 status = -7
256       RETURN
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
270   260 bound = 1.0D0
271   270 status = 3
272       RETURN
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
286   300 bound = 1.0D0
287   310 status = 4
288       RETURN
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
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
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
319   350     fx = ccum - q
320   360     CALL dinvr(status,s,fx,qleft,qhi)
321           GO TO 340
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
329   380     status = 2
330           bound = xn
331   390     CONTINUE
332   400     CONTINUE
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
348   420     fx = ccum - q
349   430     CALL dinvr(status,xn,fx,qleft,qhi)
350           GO TO 410
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
358   450     status = 2
359           bound = inf
360   460     CONTINUE
361   470     CONTINUE
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
379   490     GO TO 530
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
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
398   540     status = 2
399           bound = 1.0D0
400   550     CONTINUE
401   560 END IF
403       RETURN
405       END