!******************************************************** !* Calculate Incomplete Beta Function Ix(a,b) * !* ---------------------------------------------------- * !* SAMPLE RUN: * !* * !* Incomplete Beta Function Ix(a,b) * !* * !* A = 0.500000000000000 * !* B = 5.00000000000000 * !* x = 0.200000000000000 * !* * !* Y = 0.855072399728097 * !* * !* ---------------------------------------------------- * !* Reference: * !* "Numerical Recipes, By W.H. Press, B.P. Flannery, * !* S.A. Teukolsky and T. Vetterling, Cambridge * !* University Press, 1986" [BIBLI 08]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !******************************************************** PROGRAM TEST_BETAI real*8 a,b,x,y, BETAI a=0.5d0; b=5.d0; x=2.d-01 y = BETAI(a,b,x) print *,' ' print *,' Incomplete Beta Function Ix(A,B)' print *,' ' print *,' A =', a print *,' B =', b print *,' x =', x print *,' ' print *,' Y =', y print *,' ' END Real*8 Function BETAI(A,B,X) ! returns the incompleteBeta function Ix(a,b) real*8 A,B,X, BT, BETACF, GAMMLN IF(X.LT.0.d0.OR.X.GT.1.d0) Then print *,' BETAI: Bad argument X (must be 0<=X<=1).' RETURN END IF IF(X.EQ.0.d0.OR.X.EQ.1.d0) Then BT=0.d0 ELSE BT=DEXP(GAMMLN(A+B)-GAMMLN(A)-GAMMLN(B)+A*DLOG(X)+B*DLOG(1.d0-X)) END IF IF(X.LT.(A+1.d0)/(A+B+2.d0)) THEN BETAI=BT*BETACF(A,B,X)/A RETURN ELSE BETAI=1.d0-BT*BETACF(B,A,1.d0-X)/B RETURN END IF END Real*8 Function BETACF(A,B,X) ! continued fraction for incomplete Beta function, used by BETAI integer, parameter :: ITMAX=100 real*8, parameter :: EPS=3.d-7 real*8 A,B,X, AM,BM,AZ,QAB,QAP,QAM,BZ,EM,TEM,D,AP,BP,APP,BPP real*8 AOLD AM=1.d0 BM=1.d0 AZ=1.d0 QAB=A+B QAP=A+1.d0 QAM=A-1.d0 BZ=1.d0-QAB*X/QAP DO M=1, ITMAX EM=M TEM=EM+EM D=EM*(B-M)*X/((QAM+TEM)*(A+TEM)) AP=AZ+D*AM BP=BZ+D*BM D=-(A+EM)*(QAB+EM)*X/((A+TEM)*(QAP+TEM)) APP=AP+D*AZ BPP=BP+D*BZ AOLD=AZ AM=AP/BPP BM=BP/BPP AZ=APP/BPP BZ=1.d0 IF(DABS(AZ-AOLD).LT.EPS*DABS(AZ)) GOTO 1 END DO print *,' BETACF: A or B too big, or ITMAX too small.' RETURN 1 BETACF=AZ RETURN END Real*8 Function GAMMLN(XX) ! returns the value Ln(Gamma(XX) for XX>0. real*8 XX,COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF,STP /76.18009173d0,-86.50532033d0,24.01409822d0, & -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/ DATA HALF,ONE,FPF /0.5d0,1.d0,5.5d0/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*DLOG(TMP)-TMP SER=ONE DO J=1, 6 X=X+ONE SER=SER+COF(J)/X END DO GAMMLN=TMP+DLOG(STP*SER) RETURN END ! end of file ibeta.f90