! Last change: DOS 29 Jul 2000 11:45 am ! *** copyright 2000 *** ! *** filename bratio.f95 *** John F. Monahan ** ! ********************** program pbratio ! test problems for computing the incomplete beta ratio ! implicit none real a,b,x,beta,bratio integer j1,j2,j3 ! 21 format(5f12.6) 22 format(' a,b,x',3f6.2,6x,'beta',e20.8) ! open( unit=6, file='bratio.out' ) ! ! get values of a, b, x over a wide range do j1 = 1,5 a = .5 + (j1-1)*j1*2.5 do j2 = 1,5 b = .6 + (j2-1)*j2*2.5 do j3 = 1,10 x = .05 + .1*(j3-1) beta = bratio(x,a,b) write(6,22) a,b,x,beta ! done end do ! loop on j3 end do ! loop on j2 end do ! loop on j1 stop end program pbratio real function bratio(x,a,b) ! computes incomplete beta ratio ! defined as the integral from 0 to x of ! (a-1) (b-1) ! gamma(a+b) * x * (1-x) dx ! -------------------------------------- ! gamma(a) * gamma(b) ! ! for values ! a > 0, b > 0, and 0 < x < 1 ! ! if x < 0 then bratio = 0 ! if x > 1 then bratio = 1 ! if a or b < 0 then bratio = -1 *** error *** ! if convergence failure then bratio = -2 *** problem *** ! ! ! two main methods ! 1) power series ! 2) continued fraction due to DiDonato & Morris (1992) ! ACM TOMS, Volume 18, pp. 360 ff ! ! also employs ! 1) inversion/reflection formula ! bratio(x,a,b) = 1 - bratio(1-x,b,a) ! 2) reduction formula ! bratio(x,a,b) = ugly expression + bratio(x,a+1,b-1) ! ! J F Monahan (June 1998) Dept of Statistics, NC State U, Raleigh USA ! Recoded February, April 2000 for Fortran 95 implicit none real x,a,b,xc,yc,ac,bc,comp,term,w,s,brlog,rg,xrat,rlam, & & fj,alphj,betaj,eps real algama integer invert,k,j ! convergence criterion for power series method data eps/ 1.e-8/ ! check parameter values bratio = -1. if( (a .le. 0) .or. (b .le. 0) ) return bratio = 1. if( x .gt. 1. ) return bratio = 0. if( x .lt. 0. ) return ! ! ! main strategy ! ! both a,b > 20 then compare x to a/(a+b) ! if x > a/(a+b) then invert ! then use continued fraction method ! both a,b < 20 then compare x to 1/2 ! if x > 1/2 then invert ! then increase a, reduce b until b < 1 ! then use power series method ! one < 20, one > 20 then compare x to average of 1/2 and a/(a+b) ! if x > ave then invert ! then increase a, reduce b at most 20 times ! until both newa, newb > 20 then use continued fraction ! or until b < 1 then use power series method ! ! determine whether to invert (use reflection) ! first what to compare to comp = ( a / (a+b) + .5 ) / 2. if( (a .lt. 20.) .and. (b .lt. 20.) ) comp = .5 if( (a .gt. 20.) .and. (b .gt. 20.) ) comp = a / (a+b) ! save if not inverting invert = 0 ac = a bc = b xc = x yc = 1. - x if( x .ge. comp ) then invert = 1 ac = b bc = a xc = 1. - x yc = x end if ! ( x .ge. comp ) ! s = 0. rg = 0. ! ! set up coefficients if use reduction formula if( bc .gt. 1. ) then brlog = ac*log(xc) + (bc-1.)*log(1.-xc) & & + algama(ac+bc) - algama(ac+1.) - algama(bc) if( brlog .gt. -50. ) rg = exp(brlog) ! ! set up for raising a and dropping b term = 1. xrat = xc / yc end if ! ( bc .gt. 1. ) ! do we consider reduction or power series? ! if( (ac .le. 20.) .or. (bc .le. 20.) ) then ! ! increase a, reduce b step until (bc .gt. 1.) ! do while ( bc .gt. 1. ) s = s + term ! increase a and drop b ac = ac + 1. bc = bc - 1. term = term * xrat * bc / ac end do ! while( bc .gt. 1. ) ! bratio = s * rg end if ! ( (ac .le. 20.) .or. (bc .le. 20.) ) ! if b is an integer -- finish up exactly if ( bc .eq. 1. ) then bratio = bratio + xc**ac ! subtract if we've inverted if( invert .eq. 1 ) bratio = 1. - bratio return end if ! ( bc .eq. 1. ) if ( bc .lt. 1. ) then ! ! now use power series method ! real b is less than one ! and x is relatively small ! beta function brlog = ac*log(xc) + algama(ac+bc) - algama(ac) - algama(bc) rg = 0. if( brlog .gt. -50 ) rg = exp(brlog) ! s = 1. w = ac do k = 1,20 w = (real(k)-bc) * xc * w / real(k) term = w / (ac + real(k)) s = s + term if( abs(term) .lt. eps*s ) exit end do if( k .gt. 20 ) then ! convergence failure bratio = -2. return end if ! ( k .gt. 20 ) ! finish power series method bratio = bratio + s * rg / ac ! subtract if we've inverted if( invert .eq. 1 ) bratio = 1. - bratio return end if ! ( bc .lt. 1. ) ! continued fraction method ! if ( bc .gt. 1. ) then ! beta function brlog = ac*log(xc) + bc*log(yc) + algama(ac+bc) & & - algama(ac) - algama(bc) rg = 0. if( brlog .gt. -50 ) rg = exp(brlog) ! constant rlam = ac - (ac+bc)*xc ! compute the continued fraction s = 0. do j = 1,40 ! reverse order fj = real(40 + 1 - j) ! first the numerator alphj = (ac+fj-1.)*(ac+bc+fj-1.)*fj*(bc-fj)*xc*xc & & / ( (ac+fj+fj-1.)*(ac+fj+fj-1.) ) ! ! then the denominators betaj = fj + fj*(bc-fj)*xc/(ac+fj+fj-1.) & & + (ac+fj)*(rlam+1.+fj*(1.+yc)) / (ac+fj+fj+1.) ! s = alphj / ( betaj + s ) end do ! loop on j ! s was ratio above, now denom s = ac*(rlam + 1.) / (ac + 1.) + s bratio = bratio + ( rg / s ) ! subtract if we've inverted if( invert .eq. 1 ) bratio = 1. - bratio return end if ! ( bc .gt. 1. ) end function bratio REAL FUNCTION ALGAMA(X) ! COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS ! X .LE. 0. RETURNS -1.0 IMPLICIT NONE REAL, INTENT(IN) :: X REAL, DIMENSION(5) :: P,Q ! COEFS FOR RATIONAL FUNCTION REAL, DIMENSION(3) :: R ! COEFS FOR STIRLING CORRECTION REAL XN,A,V REAL, PARAMETER :: LSQR2PI = .918938533 DATA P/ -.510046056E1, -.368309363E2, -.664664448E2, & & -.195692802E3, -.185072903E1 / DATA Q/ -.114321842E2, .368848969E2, .16269403E2, & & -.195692802E3, 1.0 / DATA R/ -.277765545E-2, .833333330E-1, .77783067E-3 / ! USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS' ! ED. BY JF HART, SIAM SERIES, WILEY(1968) ! 27 MAR 74 JFM W/RWS,RM ! ** MODIFICATIONS JUNE 1988, 1993 ** ! RECODED OCTOBER 1999 FOR FORTRAN 95 IF( X .LE. 0. ) THEN ALGAMA = -1.0 ! NOT DESIGNED FOR NEGATIVE ARGUMENTS RETURN ENDIF ! FOR LARGE X (>8) USE STIRLING ROUTE IF( X .GT. 8. ) THEN ! LARGE X -- USE STIRLING WITH CORRECTION XN = X * X ! CORRECTION FOR STIRLING V = R(2) + ( R(3)/XN + R(1) )/XN ALGAMA = (X-0.5)*LOG(X) - X + V/X + LSQR2PI ELSE ! FOR SMALL X, USE REPRODUCTION FORMULA ! TO GET NUMBER BETWEEN 2 AND 3 V = 1.0 XN = X DO WHILE( XN .LT. 2. ) V = V/XN XN = XN + 1. END DO ! WHILE( XN .LT. 2. ) DO WHILE( XN .GT. 3. ) XN = XN - 1. V = V * XN END DO ! WHILE( XN .GT. 3. ) A = XN-2.0 ! RATIONAL FUNCTION APPROXIMATION V = V*(P(4)+A*(P(3)+A*(P(2)+A*(P(1)+A*P(5))))) / & & (Q(4)+A*(Q(3)+A*(Q(2)+A*(Q(1)+A*Q(5))))) ALGAMA = LOG(V) ENDIF RETURN END FUNCTION ALGAMA ! *** end of filename bratio.f95 *********************