FUNCTION t_quantile(p, ndf, normdev) RESULT(fn_val) ! Calculates the two-tail quantiles of Student's t-distribution. ! The user must supply a function `normdev' to return quantiles of the ! the normal distribution, such as Applied Statistics algorithm AS 241. ! Translation from Algol by Alan Miller, CSIRO Division of Mathematics ! & Statistics, Clayton, Victoria 3169, Australia of: ! Algorithm 396: Student's t-quantiles by G.W. Hill ! Comm. A.C.M., vol.13(10), 619-620, October 1970 ! N.B. The number of degrees of freedom is not necessarily an integer ! Arguments: ! P double precision Area under both tails ! ndf double precision Number of degrees of freedom ! normal double precision External function to return quantiles ! of 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) :: p, ndf REAL (dp) :: fn_val INTERFACE FUNCTION normdev(p) RESULT(x) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p REAL (dp) :: x END FUNCTION normdev END INTERFACE ! Local variables REAL (dp) :: a, b, c, d, prob, x, y REAL (dp), PARAMETER :: half_pi = 1.5707963267949D0, eps = 1.d-12, one = 1.d0 IF (ndf < one .OR. p > one .OR. p <= 0.d0) THEN WRITE(*, *) 'Error in t_quantile: ndf < 1 or p not between 0 and 1' RETURN END IF IF ( ABS(ndf - 2.d0) < eps ) THEN fn_val = SQRT(2.d0 / ( p * (2.d0 - p) ) - 2.d0) RETURN ELSE IF (ndf < one+eps) THEN prob = p * half_pi fn_val = COS(prob) / SIN(prob) RETURN ELSE a = one / (ndf - 0.5D0) b = 48.d0 / a**2 c = ((20700.d0 * a / b - 98.d0) * a - 16.d0) * a + 96.36D0 d = ((94.5D0 / (b + c) - 3.d0) / b + one) * SQRT(a * half_pi)* ndf x = d * p y = x ** (2.d0 / ndf) IF (y > 0.05D0 + a) THEN ! Asymptotic inverse expansion about normal x = normdev(0.5 * p) y = x**2 IF (ndf < 5.d0) c = c + 0.3D0 * (ndf - 4.5D0) * (x + 0.6D0) c = (((0.05D0 * d * x - 5.d0) * x - 7.d0) * x - 2.d0) * x + b+ c y = (((((0.4D0*y + 6.3D0) * y + 36.d0) * y + 94.5D0) / c - y & - 3.d0) / b + one) * x y = a * y**2 IF (y > 0.002D0) THEN y = EXP(y) - one ELSE y = 0.5D0 * y**2 + y END IF ELSE y = ((one / (((ndf + 6.d0) / (ndf * y) - 0.089D0 * d - & 0.822D0) * (ndf + 2.d0) * 3.d0) + 0.5D0 / (ndf + 4.d0)) & * y - one) * (ndf + one) / (ndf + 2.d0) + one / y END IF END IF fn_val = SQRT(ndf * y) RETURN END FUNCTION t_quantile