SUBROUTINE rhohat(rho, y, tol) !********************************************************************* ! Find RHOHAT (ML estimator) of gamma density by Y ! Requires subroutine POLG(X, PSI, POL, N, ICODE) ! INPUT: Y = LN(XBAR / Geometric mean) = LN(RHO) - PSI(RHO) ! TOL = desired accuracy ! OUTPUT: RHO = ML estimate of RHO ! N.B. The scale parameter for the gamma density is then obtained by ! equating the population (arithmetic) mean to XBAR = sample mean. ! This version, which is compatible with Lahey's ELF90 compiler, ! is by Alan Miller: amiller @ bigpond.net.au ! URL: http://users.bigpond.net.au/amiller ! Latest revision - 27 August 1997 !********************************************************************* IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: y, tol REAL (dp), INTENT(OUT) :: rho ! Local variables REAL (dp) :: r1, di REAL (dp), PARAMETER :: one = 1._dp, four = 4._dp, three = 3._dp INTEGER :: i INTERFACE SUBROUTINE polg(x, psi, pol, n, icode) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: x INTEGER, INTENT(IN), OPTIONAL :: n REAL (dp), INTENT(OUT) :: psi REAL (dp), INTENT(OUT), OPTIONAL :: pol(:) INTEGER, INTENT(OUT), OPTIONAL :: icode END SUBROUTINE polg END INTERFACE IF (y <= EPSILON(one)) THEN WRITE(*, *) 'Error: argument y to routine RHOHAT must be > 0' RETURN END IF ! Initial RHO is Thom's estimator rho = (one + SQRT(one + four * y / three)) / (four * y) DO i = 1, 100 r1 = rho CALL polg(r1, di) rho = r1 * (LOG(r1) - di) / y IF (ABS(r1 - rho) < tol) EXIT END DO RETURN END SUBROUTINE rhohat SUBROUTINE polg(x, psi, pol, n, icode) ! This program computes PSI(X) & its first N derivatives. !********************************************************************* ! Input parameter X must be positive and < 4000 ! N integer - highest derivative required >= 0 ! Output parameter PSI = 1st derivative of log (gamma (X) ) ! POL dimension at least N contains the derivatives ! ICODE = 0 for normal exit ! = 1 if X <= 0 ! = 2 if N < 0 or N > 20 !********************************************************************* ! BER(S) = the Bernouilli number of order 2S (apart from signs) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: x INTEGER, INTENT(IN), OPTIONAL :: n REAL (dp), INTENT(OUT) :: psi REAL (dp), INTENT(OUT), OPTIONAL :: pol(:) INTEGER, INTENT(OUT), OPTIONAL :: icode ! Local variables REAL (dp) :: x2, w(8), di, dj, dum, ds, sgn REAL (dp), PARAMETER :: ber(8) = (/ 0.16666666666666667D+00, & 0.33333333333333333D-01, 0.23809523809523810D-01, 0.33333333333333333D-01, & 0.75757575757575758D-01, 0.25311355311355311D+00, 0.11666666666666667D+01, & 0.70921568627450980D+01 /) ! Values above are 1/6, 1/30, 1/42, 1/30, 5/66, 691/2730, 7/6 & 3617/510. INTEGER :: i, ic, in, indx, j, jj, nn ! Check for valid input parameters X and N. ic = 0 IF (x < 0.01D0 .OR. x > 4000.D0) ic = 1 IF (PRESENT(n)) THEN IF (n < 0 .OR. n > 20) ic = 2 nn = n ELSE nn = 0 END IF IF (PRESENT(icode)) icode = ic IF (ic > 0) RETURN ! Increment X to a value > 40 indx = 0 x2 = x DO i = 1, 40 IF (x2 > 40.D0) EXIT indx = indx + 1 x2 = x2 + 1.D0 END DO IF (x2 > 500.D0) THEN in = 3 ELSE IF (x2 > 200.D0) THEN in = 4 ELSE IF (x2 > 100.D0) THEN in = 6 ELSE in = 8 END IF ! Calculate polygamma function using table 6.4.11 on page 260 of ! Bowman, K.O. & Shenton, L.R. (1988), Properties of Estimators for ! the Gamma Distribution, Marcel Dekker. ! X2 = 40 + fractional part of X IF (nn == 0) GO TO 40 dum = 1.D0 / x2 sgn = 1.D0 DO i = 1, nn di = i pol(i) = dum * (1.D0 + di * 0.5D0 / x2) dum = - dum * di / x2 DO j = 1, in dj = j IF (i == 1) THEN w(j) = sgn * ber(j) / x2 ** (2*j+1) sgn = - sgn ELSE w(j) = w(j) * (1.D0 - di - 2.D0 * dj) / x2 END IF pol(i) = pol(i) + w(j) END DO END DO 40 psi = LOG(x2) - 0.5D0 / x2 sgn = - 1.D0 DO i = 1, in di = 2 * i dum = ber(i) / x2 ** (2*i) / di psi = psi + sgn * dum sgn = - sgn END DO ! Calculate desirable polygamma function using 6.4.6 from page 260. IF (indx == 0) RETURN IF (PRESENT(pol)) THEN sgn = 1.D0 DO i = 1, nn ds = 0.D0 DO j = 1, indx dj = indx - j dum = 1.D0 / (x + dj) dj = dum DO jj = 1, i di = jj dj = dj * di * dum END DO ds = ds + dj END DO pol(i) = pol(i) + sgn * ds sgn = - sgn END DO END IF DO i = 1, indx dum = i - 1 psi = psi - 1.D0 / (x + dum) END DO RETURN END SUBROUTINE polg PROGRAM drive IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp) :: y, rho, tol = 1.D-04 INTERFACE SUBROUTINE rhohat(rho, y, tol) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: y, tol REAL (dp), INTENT(OUT) :: rho END SUBROUTINE rhohat END INTERFACE DO WRITE(*, '(a)', ADVANCE='NO') ' Enter value of y: ' READ(*, *) y CALL rhohat(rho, y, tol) WRITE(*, '(a, f10.6)') ' RHO = ', rho WRITE(*, *) END DO STOP END PROGRAM drive