FUNCTION student(t, ndf, normal) RESULT(fn_val) ! Calculates the two-tailed probability of Student's t. ! The user must supply a function `normal' to calculate the area under the ! normal distribution curve, such as Applied Statistics algorithm AS 66. ! Translation from Algol by Alan Miller, CSIRO Division of Mathematics ! & Statistics, Clayton, Victoria 3169, Australia of: ! Algorithm 395: Student's t-distribution by G.W. Hill ! Comm. A.C.M., vol.13(10), 617-619, October 1970 ! N.B. The number of degrees of freedom is not necessarily an integer ! Arguments: ! t double precision Value of Student's t ! ndf double precision Number of degrees of freedom ! normal double precision External function to return the ! area under the standard normal curve ! ier integer Error indicator = 0 for normal exit ! = 1 if ndf < 1 ! ELF90-compatible version by Alan Miller ! Latest revision - 7 August 1997 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: t, ndf REAL (dp) :: fn_val INTERFACE FUNCTION normal(x) RESULT(fx) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: x REAL (dp) :: fx END FUNCTION normal END INTERFACE ! Local variables REAL (dp) :: a, b, n, t2, y, z REAL (dp), PARAMETER :: eps = 1.d-12, one = 1.d0, twoonpi = 0.63661977236758D0 INTEGER :: j IF (ndf < one) THEN WRITE(*, *) 'Error in student: Number of degrees of freedom < 1' RETURN END IF t2 = t*t y = t2/ndf b = one + y z = one IF ( (ndf - INT(ndf)) > eps .OR. (ndf > 20.d0 - eps .AND. ndf > t2) & .OR. ndf > 200.d0 ) THEN ! Asymptotic series for large or non-integer ndf IF (y > 1.d-06) y = LOG(b) a = ndf - 0.5D0 b = 48.d0 * a**2 y = a * y y = (((((-0.4D0*y - 3.3D0)*y - 24.d0)*y - 85.5D0) / & (0.8D0*y**2 + 100.d0 + b) + y +3.d0) / b + one) * SQRT(y) fn_val = 2.d0 * normal(-y) RETURN ELSE IF ( ndf < 20.d0 .AND. t2 < 4.d0) THEN ! Nested summation of cosine series y = SQRT(y) a = y IF (ndf == one) a = 0.d0 n = ndf ! Tail series expansion for large t-values ELSE a = SQRT(b) y = a * ndf j = 0 10 j = j + 2 IF (ABS(a-z) > eps) THEN z = a y = y * (j - 1) / (b * j) a = a + y / (ndf + j) GO TO 10 END IF n = ndf + 2.d0 z = 0.d0 y = 0.d0 a = - a END IF ! This is the `loop' from the Algol code. ! It required jumping between different parts of an IF () THEN ! ELSE IF () .. block. 20 n = n - 2.d0 IF (n > 1.d0) THEN a = (n - 1.d0) * a / (b * n) + y GO TO 20 END IF IF ( ABS(n) < eps ) THEN a = a / SQRT(b) ELSE a = (ATAN(y) + a/b) * twoonpi END IF fn_val = z - a RETURN END FUNCTION student