C C .................................................................. C C SUBROUTINE BDTR C C PURPOSE C COMPUTES P(X) = PROBABILITY THAT THE RANDOM VARIABLE U, C DISTRIBUTED ACCORDING TO THE BETA DISTRIBUTION WITH C PARAMETERS A AND B, IS LESS THAN OR EQUAL TO X. F(A,B,X), C THE ORDINATE OF THE BETA DENSITY AT X, IS ALSO COMPUTED. C C USAGE C CALL BDTR(X,A,B,P,D,IER) C C DESCRIPTION OF PARAMETERS C X - INPUT SCALAR FOR WHICH P(X) IS COMPUTED. C A - BETA DISTRIBUTION PARAMETER (CONTINUOUS). C B - BETA DISTRIBUTION PARAMETER (CONTINUOUS). C P - OUTPUT PROBABILITY. C D - OUTPUT DENSITY. C IER - RESULTANT ERROR CODE WHERE C IER= 0 --- NO ERROR C IER=-1,+1 CDTR HAS BEEN CALLED AND AN ERROR HAS C OCCURRED. SEE CDTR. C IER=-2 --- AN INPUT PARAMETER IS INVALID. X IS LESS C THAN 0.0 OR GREATER THAN 1.0, OR EITHER A OR C B IS LESS THAN 0.5 OR GREATER THAN 10**(+5). C P AND D ARE SET TO -1.7E38. 0 C IER=+2 --- INVALID OUTPUT. P IS LESS THAN ZERO OR C GREATER THAN ONE. P IS SET TO 1.7E38. 0 C C REMARKS C SEE MATHEMATICAL DESCRIPTION. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C DLGAM C NDTR C CDTR C C METHOD C REFER TO R.E. BARGMANN AND S.P. GHOSH, STATISTICAL C DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE, C IBM RESEARCH REPORT RC-1094, 1963. C C .................................................................. C SUBROUTINE BDTR(X,A,B,P,D,IER) DOUBLE PRECISION XX,DLXX,DL1X,AA,BB,G1,G2,G3,G4,DD,PP,XO,FF,FN, 1XI,SS,CC,RR,DLBETA C C TEST FOR VALID INPUT DATA C IF(A-(.5-1.E-5)) 640,10,10 10 IF(B-(.5-1.E-5)) 640,20,20 20 IF(A-1.E+5) 30,30,640 30 IF(B-1.E+5) 40,40,640 40 IF(X) 640,50,50 50 IF(1.-X) 640,60,60 C C COMPUTE LOG(BETA(A,B)) C 60 AA=DBLE(A) BB=DBLE(B) CALL DLGAM(AA,G1,IOK) CALL DLGAM(BB,G2,IOK) CALL DLGAM(AA+BB,G3,IOK) DLBETA=G1+G2-G3 C C TEST FOR X NEAR 0.0 OR 1.0 C IF(X-1.E-8) 80,80,70 70 IF((1.-X)-1.E-8) 130,130,140 80 P=0.0 IF(A-1.) 90,100,120 90 D=1.7E38 GO TO 660 100 DD=-DLBETA IF(DD+1.68D02) 120,120,110 110 DD=DEXP(DD) D=SNGL(DD) GO TO 660 120 D=0.0 GO TO 660 130 P=1.0 IF(B-1.) 90,100,120 C C SET PROGRAM PARAMETERS C 140 XX=DBLE(X) DLXX=DLOG(XX) DL1X=DLOG(1.D0-XX) XO=XX/(1.D0-XX) ID=0 C C COMPUTE ORDINATE C DD=(AA-1.D0)*DLXX+(BB-1.D0)*DL1X-DLBETA IF(DD-1.68D02) 150,150,160 150 IF(DD+1.68D02) 170,170,180 160 D=1.7E38 0 GO TO 190 170 D=0.0 GO TO 190 180 DD=DEXP(DD) D=SNGL(DD) C C A OR B OR BOTH WITHIN 1.E-8 OF 1.0 C 190 IF(ABS(A-1.)-1.E-8) 200,200,210 200 IF(ABS(B-1.)-1.E-8) 220,220,230 210 IF(ABS(B-1.)-1.E-8) 260,260,290 220 P=X GO TO 660 230 PP=BB*DL1X IF(PP+1.68D02) 240,240,250 240 P=1.0 GO TO 660 250 PP=DEXP(PP) PP=1.D0-PP P=SNGL(PP) GO TO 600 260 PP=AA*DLXX IF(PP+1.68D02) 270,270,280 270 P=0.0 GO TO 660 280 PP=DEXP(PP) P=SNGL(PP) GO TO 600 C C TEST FOR A OR B GREATER THAN 1000.0 C 290 IF(A-1000.) 300,300,310 300 IF(B-1000.) 330,330,320 310 XX=2.D0*AA/XO XS=SNGL(XX) AA=2.D0*BB DF=SNGL(AA) CALL CDTR(XS,DF,P,DUMMY,IER) P=1.0-P GO TO 670 320 XX=2.D0*BB*XO XS=SNGL(XX) AA=2.D0*AA DF=SNGL(AA) CALL CDTR(XS,DF,P,DUMMY,IER) GO TO 670 C C SELECT PARAMETERS FOR CONTINUED FRACTION COMPUTATION C 330 IF(X-.5) 340,340,380 340 IF(AA-1.D0) 350,350,360 350 RR=AA+1.D0 GO TO 370 360 RR=AA 370 DD=DLXX/5.D0 DD=DEXP(DD) DD=(RR-1.D0)-(RR+BB-1.D0)*XX*DD +2.D0 IF(DD) 420,420,430 380 IF(BB-1.D0) 390,390,400 390 RR=BB+1.D0 GO TO 410 400 RR=BB 410 DD=DL1X/5.D0 DD=DEXP(DD) DD=(RR-1.D0)-(AA+RR-1.D0)*(1.D0-XX)*DD +2.D0 IF(DD) 430,430,420 420 ID=1 FF=DL1X DL1X=DLXX DLXX=FF XO=1.D0/XO FF=AA AA=BB BB=FF G2=G1 C C TEST FOR A LESS THAN 1.0 C 430 FF=0.D0 IF(AA-1.D0) 440,440,470 440 CALL DLGAM(AA+1.D0,G4,IOK) DD=AA*DLXX+BB*DL1X+G3-G2-G4 IF(DD+1.68D02) 460,460,450 450 FF=FF+DEXP(DD) 460 AA=AA+1.D0 C C COMPUTE P USING CONTINUED FRACTION EXPANSION C 470 FN=AA+BB-1.D0 RR=AA-1.D0 II=80 XI=DFLOAT(II) SS=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI)) SS=SS*XO DO 480 I=1,79 II=80-I XI=DFLOAT(II) DD=(XI*(FN+XI))/((RR+2.D0*XI+1.D0)*(RR+2.D0*XI)) DD=DD*XO CC=((BB-XI)*(RR+XI))/((RR+2.D0*XI-1.D0)*(RR+2.D0*XI)) CC=CC*XO SS=CC/(1.D0+DD/(1.D0-SS)) 480 CONTINUE SS=1.D0/(1.D0-SS) IF(SS) 650,650,490 490 CALL DLGAM(AA+BB,G1,IOK) CALL DLGAM(AA+1.D0,G4,IOK) CC=G1-G2-G4+AA*DLXX+(BB-1.D0)*DL1X PP=CC+DLOG(SS) IF(PP+1.68D02) 500,500,510 500 PP=FF GO TO 520 510 PP=DEXP(PP)+FF 520 IF(ID) 540,540,530 530 PP=1.D0-PP 540 P=SNGL(PP) C C SET ERROR INDICATOR C IF(P) 550,570,570 550 IF(ABS(P)-1.E-7) 560,560,650 560 P=0.0 GO TO 660 570 IF(1.-P) 580,600,600 580 IF(ABS(1.-P)-1.E-7) 590,590,650 590 P=1.0 GO TO 660 600 IF(P-1.E-8) 610,610,620 610 P=0.0 GO TO 660 620 IF((1.0-P)-1.E-8) 630,630,660 630 P=1.0 GO TO 660 640 IER=-2 D=-1.7E38 0 P=-1.7E38 0 GO TO 670 650 IER=+2 P= 1.7E38 0 GO TO 670 660 IER=0 670 RETURN END