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