FUNCTION ppchi2(p, v, g) RESULT(fn_val) ! N.B. Argument IFAULT has been removed. ! This version by Alan Miller ! amiller @ bigpond.net.au ! Latest revision - 27 October 2000 ! Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35 ! To evaluate the percentage points of the chi-squared ! probability distribution function. ! p must lie in the range 0.000002 to 0.999998, ! v must be positive, ! g must be supplied and should be equal to ln(gamma(v/2.0)) ! Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991) ! which should eliminate the need for the limited range for p above, ! though these limits have not been removed from the routine. ! If IFAULT = 4 is returned, the result is probably as accurate as ! the machine will allow. ! Auxiliary routines required: PPND = AS 111 (or AS 241) and GAMMAD = AS 239. IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p REAL (dp), INTENT(IN) :: v REAL (dp), INTENT(IN) :: g REAL (dp) :: fn_val INTERFACE FUNCTION gammad(x, p) RESULT(fn_val) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: x, p REAL (dp) :: fn_val END FUNCTION gammad SUBROUTINE ppnd16 (p, normal_dev, ifault) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: p INTEGER, INTENT(OUT) :: ifault REAL (dp), INTENT(OUT) :: normal_dev END SUBROUTINE ppnd16 END INTERFACE ! Local variables REAL (dp) :: a, b, c, p1, p2, q, s1, s2, s3, s4, s5, s6, t, x, xx INTEGER :: i, if1 INTEGER, PARAMETER :: maxit = 20 REAL (dp), PARAMETER :: aa = 0.6931471806_dp, e = 0.5e-06_dp, & pmin = 0.000002_dp, pmax = 0.999998_dp, & zero = 0.0_dp, half = 0.5_dp, one = 1.0_dp, & two = 2.0_dp, three = 3.0_dp, six = 6.0_dp, & c1 = 0.01_dp, c2 = 0.222222_dp, c3 = 0.32_dp, & c4 = 0.4_dp, c5 = 1.24_dp, c6 = 2.2_dp, & c7 = 4.67_dp, c8 = 6.66_dp, c9 = 6.73_dp, & c10 = 13.32_dp, c11 = 60.0_dp, c12 = 70.0_dp, & c13 = 84.0_dp, c14 = 105.0_dp, c15 = 120.0_dp, & c16 = 127.0_dp, c17 = 140.0_dp, c18 = 175.0_dp, & c19 = 210.0_dp, c20 = 252.0_dp, c21 = 264.0_dp, & c22 = 294.0_dp, c23 = 346.0_dp, c24 = 420.0_dp, & c25 = 462.0_dp, c26 = 606.0_dp, c27 = 672.0_dp, & c28 = 707.0_dp, c29 = 735.0_dp, c30 = 889.0_dp, & c31 = 932.0_dp, c32 = 966.0_dp, c33 = 1141.0_dp, & c34 = 1182.0_dp, c35 = 1278.0_dp, c36 = 1740.0_dp, & c37 = 2520.0_dp, c38 = 5040.0_dp ! Test arguments and initialise fn_val = -one IF (p < pmin .OR. p > pmax) THEN WRITE(*, *) 'Error in PPCHI2: p must be between 0.000002 & 0.999998' RETURN END IF IF (v <= zero) THEN WRITE(*, *) 'Error in PPCHI2: Number of deg. of freedom <= 0' RETURN END IF xx = half * v c = xx - one ! Starting approximation for small chi-squared IF (v < -c5 * LOG(p)) THEN fn_val = (p * xx * EXP(g + xx * aa)) ** (one/xx) IF (fn_val < e) GO TO 6 GO TO 4 END IF ! Starting approximation for v less than or equal to 0.32 IF (v > c3) GO TO 3 fn_val = c4 a = LOG(one-p) 2 q = fn_val p1 = one + fn_val * (c7+fn_val) p2 = fn_val * (c9 + fn_val * (c8 + fn_val)) t = -half + (c7 + two * fn_val) / p1 - (c9 + fn_val * (c10 + three * fn_val)) / p2 fn_val = fn_val - (one - EXP(a + g + half * fn_val + c * aa) * p2 / p1) / t IF (ABS(q / fn_val - one) > c1) GO TO 2 GO TO 4 ! Call to algorithm AS 241 - note that p has been tested above. 3 CALL ppnd16(p, x, if1) ! Starting approximation using Wilson and Hilferty estimate p1 = c2 / v fn_val = v * (x * SQRT(p1) + one - p1) ** 3 ! Starting approximation for p tending to 1 IF (fn_val > c6 * v + six) fn_val = -two * (LOG(one-p) - c * LOG(half * fn_val) + g) ! Call to algorithm AS 239 and calculation of seven term Taylor series 4 DO i = 1, maxit q = fn_val p1 = half * fn_val p2 = p - gammad(p1, xx) t = p2 * EXP(xx * aa + g + p1 - c * LOG(fn_val)) b = t / fn_val a = half * t - b * c s1 = (c19 + a * (c17 + a * (c14 + a * (c13 + a * (c12 + c11 * a))))) / c24 s2 = (c24 + a * (c29 + a * (c32 + a * (c33 + c35 * a)))) / c37 s3 = (c19 + a * (c25 + a * (c28 + c31 * a))) / c37 s4 = (c20 + a * (c27 + c34 * a) + c * (c22 + a * (c30 + c36 * a))) / c38 s5 = (c13 + c21 * a + c * (c18 + c26 * a)) / c37 s6 = (c15 + c * (c23 + c16 * c)) / c38 fn_val = fn_val + t * (one + half * t * s1 - b * c * (s1 - b * & (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6)))))) IF (ABS(q / fn_val - one) > e) RETURN END DO WRITE(*, *) 'Error in PPCHI2: Max. number of iterations exceeded' 6 RETURN END FUNCTION ppchi2