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