1       SUBROUTINE cdfbet(which,p,q,x,y,a,b,status,bound)
2 C**********************************************************************
3 C
4 C      SUBROUTINE CDFBET( WHICH, P, Q, X, Y, A, B, STATUS, BOUND )
5 C               Cumulative Distribution Function
6 C                         BETa Distribution
7 C
8 C
9 C                              Function
10 C
11 C
12 C     Calculates any one parameter of the beta distribution given
13 C     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 X,Y,A and B
23 C               iwhich = 2 : Calculate X and Y from P,Q,A and B
24 C               iwhich = 3 : Calculate A from P,Q,X,Y and B
25 C               iwhich = 4 : Calculate B from P,Q,X,Y and A
26 C
27 C                    INTEGER WHICH
28 C
29 C     P <--> The integral from 0 to X of the chi-square
30 C            distribution.
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     X <--> Upper limit of integration of beta density.
40 C            Input range: [0,1].
41 C            Search range: [0,1]
42 C                    DOUBLE PRECISION X
43 C
44 C     Y <--> 1-X.
45 C            Input range: [0,1].
46 C            Search range: [0,1]
47 C            X + Y = 1.0.
48 C                    DOUBLE PRECISION Y
49 C
50 C     A <--> The first parameter of the beta density.
51 C            Input range: (0, +infinity).
52 C            Search range: [1D-300,1D300]
53 C                    DOUBLE PRECISION A
54 C
55 C     B <--> The second parameter of the beta density.
56 C            Input range: (0, +infinity).
57 C            Search range: [1D-300,1D300]
58 C                    DOUBLE PRECISION B
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 X + Y .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     Cumulative distribution function  (P)  is calculated directly by
84 C     code associated with the following reference.
85 C
86 C     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
87 C     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
88 C     Trans. Math.  Softw. 18 (1993), 360-373.
89 C
90 C     Computation of other parameters involve a seach for a value that
91 C     produces  the desired  value  of P.   The search relies  on  the
92 C     monotinicity of P with the other parameter.
93 C
94 C
95 C                              Note
96 C
97 C
98 C     The beta density is proportional to
99 C               t^(A-1) * (1-t)^(B-1)
100 C
101 C**********************************************************************
102 C     .. Parameters ..
103       DOUBLE PRECISION tol
104       PARAMETER (tol=1.0D-13)
105       DOUBLE PRECISION atol
106       PARAMETER (atol=1.0D-50)
107       DOUBLE PRECISION zero,inf
108       PARAMETER (zero=1.0D-300,inf=1.0D300)
109       DOUBLE PRECISION one
110       PARAMETER (one=1.0D0)
111 C     ..
112 C     .. Scalar Arguments ..
113       DOUBLE PRECISION a,b,bound,p,q,x,y
114       INTEGER status,which
115 C     ..
116 C     .. Local Scalars ..
117       DOUBLE PRECISION fx,xhi,xlo,cum,ccum,xy,pq
118       LOGICAL qhi,qleft,qporq
119 C     ..
120 C     .. External Functions ..
121       INTEGER vfinite
122       DOUBLE PRECISION spmpar
123       EXTERNAL spmpar
124 C     ..
125 C     .. External Subroutines ..
126       EXTERNAL cumbet,dinvr,dstinv,dstzr,dzror
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
137    10 bound = 4.0D0
138    20 status = -1
139       RETURN
141    30 IF (which.EQ.1) GO TO 70
142 C
143 C     P
144 C
145       IF (ISANAN(p).EQ.1) THEN
146          CALL RETURNANANFORTRAN(a)
147          CALL RETURNANANFORTRAN(b)
148          CALL RETURNANANFORTRAN(x)
149          CALL RETURNANANFORTRAN(y)
150          RETURN
151       ENDIF
152       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
153       IF (.NOT. (p.LT.0.0D0)) GO TO 40
154       bound = 0.0D0
155       GO TO 50
157    40 bound = 1.0D0
158    50 status = -2
159       RETURN
161    60 CONTINUE
162    70 IF (which.EQ.1) GO TO 110
163 C
164 C     Q
165 C
166       IF (ISANAN(q).EQ.1) THEN
167          CALL RETURNANANFORTRAN(a)
168          CALL RETURNANANFORTRAN(b)
169          CALL RETURNANANFORTRAN(x)
170          CALL RETURNANANFORTRAN(y)
171          RETURN
172       ENDIF
173       IF (.NOT. ((q.LT.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
174       IF (.NOT. (q.LT.0.0D0)) GO TO 80
175       bound = 0.0D0
176       GO TO 90
178    80 bound = 1.0D0
179    90 status = -3
180       RETURN
182   100 CONTINUE
183   110 IF (which.EQ.2) GO TO 150
184 C
185 C     X
186 C
187       IF (ISANAN(x).EQ.1) THEN
188          CALL RETURNANANFORTRAN(p)
189          CALL RETURNANANFORTRAN(q)
190          CALL RETURNANANFORTRAN(a)
191          CALL RETURNANANFORTRAN(b)
192          RETURN
193       ENDIF
194       IF (.NOT. ((x.LT.0.0D0).OR. (x.GT.1.0D0))) GO TO 140
195       IF (.NOT. (x.LT.0.0D0)) GO TO 120
196       bound = 0.0D0
197       GO TO 130
199   120 bound = 1.0D0
200   130 status = -4
201       RETURN
203   140 CONTINUE
204   150 IF (which.EQ.2) GO TO 190
205 C
206 C     Y
207 C
208       IF (ISANAN(y).EQ.1) THEN
209          CALL RETURNANANFORTRAN(p)
210          CALL RETURNANANFORTRAN(q)
211          CALL RETURNANANFORTRAN(a)
212          CALL RETURNANANFORTRAN(b)
213          RETURN
214       ENDIF
215       IF (.NOT. ((y.LT.0.0D0).OR. (y.GT.1.0D0))) GO TO 180
216       IF (.NOT. (y.LT.0.0D0)) GO TO 160
217       bound = 0.0D0
218       GO TO 170
220   160 bound = 1.0D0
221   170 status = -5
222       RETURN
224   180 CONTINUE
225   190 IF (which.EQ.3) GO TO 210
226 C
227 C     A
228 C
229       IF (ISANAN(a).EQ.1) THEN
230          CALL RETURNANANFORTRAN(p)
231          CALL RETURNANANFORTRAN(q)
232          CALL RETURNANANFORTRAN(x)
233          CALL RETURNANANFORTRAN(y)
234          CALL RETURNANANFORTRAN(b)
235          RETURN
236       ENDIF
237       IF (vfinite(1,a).EQ.0) a = SIGN(inf,a)
238       IF (.NOT. (a.LE.0.0D0)) GO TO 200
239       bound = 0.0D0
240       status = -6
241       RETURN
243   200 CONTINUE
244   210 IF (which.EQ.4) GO TO 230
245 C
246 C     B
247 C
248       IF (ISANAN(b).EQ.1) THEN
249          CALL RETURNANANFORTRAN(p)
250          CALL RETURNANANFORTRAN(q)
251          CALL RETURNANANFORTRAN(x)
252          CALL RETURNANANFORTRAN(y)
253          CALL RETURNANANFORTRAN(a)
254          RETURN
255       ENDIF
256       IF (vfinite(1,b).EQ.0) b = SIGN(inf,b)
257       IF (.NOT. (b.LE.0.0D0)) GO TO 220
258       bound = 0.0D0
259       status = -7
260       RETURN
262   220 CONTINUE
263   230 IF (which.EQ.1) GO TO 270
264 C
265 C     P + Q
266 C
267       pq = p + q
268       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
269      +    (3.0D0*spmpar(1)))) GO TO 260
270       IF (.NOT. (pq.LT.0.0D0)) GO TO 240
271       bound = 0.0D0
272       GO TO 250
274   240 bound = 1.0D0
275   250 status = 3
276       RETURN
278   260 CONTINUE
279   270 IF (which.EQ.2) GO TO 310
280 C
281 C     X + Y
282 C
283       xy = x + y
284       IF (.NOT. (abs(((xy)-0.5D0)-0.5D0).GT.
285      +    (3.0D0*spmpar(1)))) GO TO 300
286       IF (.NOT. (xy.LT.0.0D0)) GO TO 280
287       bound = 0.0D0
288       GO TO 290
290   280 bound = 1.0D0
291   290 status = 4
292       RETURN
294   300 CONTINUE
295   310 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
296 C
297 C     Select the minimum of P or Q
298 C
299 C
301 C
302       IF ((1).EQ. (which)) THEN
303 C
304 C     Calculating P and Q
305 C
306           CALL cumbet(x,y,a,b,p,q)
307           status = 0
309       ELSE IF ((2).EQ. (which)) THEN
310 C
311 C     Calculating X and Y
312 C
313           CALL dstzr(0.0D0,1.0D0,atol,tol)
314           IF (.NOT. (qporq)) GO TO 340
315           status = 0
316           CALL dzror(status,x,fx,xlo,xhi,qleft,qhi)
317           y = one - x
318   320     IF (.NOT. (status.EQ.1)) GO TO 330
319           CALL cumbet(x,y,a,b,cum,ccum)
320           fx = cum - p
321           CALL dzror(status,x,fx,xlo,xhi,qleft,qhi)
322           y = one - x
323 c          write(6,'(''x'',e10.3,''y='',e10.3,''sta='',i3)') x,y,status
324           GO TO 320
326   330     GO TO 370
328   340     status = 0
329           CALL dzror(status,y,fx,xlo,xhi,qleft,qhi)
330           x = one - y
331   350     IF (.NOT. (status.EQ.1)) GO TO 360
332           CALL cumbet(x,y,a,b,cum,ccum)
333           fx = ccum - q
334           CALL dzror(status,y,fx,xlo,xhi,qleft,qhi)
335           x = one - y
336           GO TO 350
338   360     CONTINUE
339   370     IF (.NOT. (status.EQ.-1)) GO TO 400
340           IF (.NOT. (qleft)) GO TO 380
341           status = 1
342           bound = 0.0D0
343           GO TO 390
345   380     status = 2
346           bound = 1.0D0
347   390     CONTINUE
348   400     CONTINUE
350       ELSE IF ((3).EQ. (which)) THEN
351 C
352 C     Computing A
353 C
354           a = 5.0D0
355           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
356           status = 0
357           CALL dinvr(status,a,fx,qleft,qhi)
358   410     IF (.NOT. (status.EQ.1)) GO TO 440
359           CALL cumbet(x,y,a,b,cum,ccum)
360           IF (.NOT. (qporq)) GO TO 420
361           fx = cum - p
362           GO TO 430
364   420     fx = ccum - q
365   430     CALL dinvr(status,a,fx,qleft,qhi)
366           GO TO 410
368   440     IF (.NOT. (status.EQ.-1)) GO TO 470
369           IF (.NOT. (qleft)) GO TO 450
370           status = 1
371           bound = zero
372           GO TO 460
374   450     status = 2
375           bound = inf
376   460     CONTINUE
377   470     CONTINUE
379       ELSE IF ((4).EQ. (which)) THEN
380 C
381 C     Computing B
382 C
383           b = 5.0D0
384           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
385           status = 0
386           CALL dinvr(status,b,fx,qleft,qhi)
387   480     IF (.NOT. (status.EQ.1)) GO TO 510
388           CALL cumbet(x,y,a,b,cum,ccum)
389           IF (.NOT. (qporq)) GO TO 490
390           fx = cum - p
391           GO TO 500
393   490     fx = ccum - q
394   500     CALL dinvr(status,b,fx,qleft,qhi)
395           GO TO 480
397   510     IF (.NOT. (status.EQ.-1)) GO TO 540
398           IF (.NOT. (qleft)) GO TO 520
399           status = 1
400           bound = zero
401           GO TO 530
403   520     status = 2
404           bound = inf
405   530     CONTINUE
406   540 END IF
408       RETURN
409 C
410       END