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