PROGRAM test_chi_squared IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) INTEGER :: ndf REAL (dp) :: chi2 INTERFACE FUNCTION chi_squared(ndf, chi2) RESULT(prob) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) INTEGER, INTENT(IN) :: ndf REAL (dp), INTENT(IN) :: chi2 REAL (dp) :: prob END FUNCTION chi_squared END INTERFACE DO WRITE(*, *) 'Enter no. of degrees of freedom: ' READ(*, *) ndf WRITE(*, *) 'Enter chi-squared value: ' READ(*, *) chi2 WRITE(*, *) 'Chi-squared probability = ', chi_squared(ndf, chi2) END DO STOP END PROGRAM test_chi_squared FUNCTION chi_squared(ndf, chi2) RESULT(prob) ! Calculate the chi-squared distribution function ! ndf = number of degrees of freedom ! chi2 = chi-squared value ! prob = probability of a chi-squared value <= chi2 (i.e. the left-hand ! tail area) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) INTEGER, INTENT(IN) :: ndf REAL (dp), INTENT(IN) :: chi2 REAL (dp) :: prob INTERFACE FUNCTION gammad(x, p) RESULT(gamma_prob) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL (dp), INTENT(IN) :: x, p REAL (dp) :: gamma_prob END FUNCTION gammad END INTERFACE ! Local variables REAL (dp) :: half = 0.5d0, x, p x = half * chi2 p = half * REAL(ndf) prob = gammad(x, p) RETURN END FUNCTION chi_squared FUNCTION gammad(x, p) RESULT(fn_val) ! ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3 ! Computation of the Incomplete Gamma Integral ! Auxiliary functions required: lngamma = logarithm of the gamma ! function, and ALNORM = algorithm AS66 ! ELF90-compatible version by Alan Miller ! amiller@bigpond.net.au ! Latest revision - 28 October 2000 ! N.B. Argument IFAULT has been removed IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), INTENT(IN) :: x, p REAL (dp) :: fn_val INTERFACE FUNCTION alnorm( x, upper ) RESULT( fn_val ) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 100) REAL(dp), INTENT(IN) :: x LOGICAL, INTENT(IN) :: upper REAL(DP) :: fn_val END FUNCTION alnorm FUNCTION lngamma(z) RESULT(lanczos) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL(dp), INTENT(IN) :: z REAL(dp) :: lanczos END FUNCTION lngamma END INTERFACE ! Local variables REAL (dp) :: pn1, pn2, pn3, pn4, pn5, pn6, arg, c, rn, a, b, an REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, two = 2.0_dp, & oflo = 1.d+37, three = 3.0_dp, nine = 9.0_dp, & tol = 1.d-14, xbig = 1.d+8, plimit = 1000.0_dp, & elimit = -88.0_dp ! EXTERNAL lngamma, alnorm fn_val = zero ! Check that we have valid values for X and P IF (p <= zero .OR. x < zero) THEN WRITE(*, *) 'AS239: Either p <= 0 or x < 0' RETURN END IF IF (x == zero) RETURN ! Use a normal approximation if P > PLIMIT IF (p > plimit) THEN pn1 = three * SQRT(p) * ((x / p) ** (one / three) + one /(nine * p) - one) fn_val = alnorm(pn1, .false.) RETURN END IF ! If X is extremely large compared to P then set fn_val = 1 IF (x > xbig) THEN fn_val = one RETURN END IF IF (x <= one .OR. x < p) THEN ! Use Pearson's series expansion. ! (Note that P is not large enough to force overflow in lngamma). ! No need to test IFAULT on exit since P > 0. arg = p * LOG(x) - x - lngamma(p + one) c = one fn_val = one a = p DO a = a + one c = c * x / a fn_val = fn_val + c IF (c <= tol) EXIT END DO arg = arg + LOG(fn_val) fn_val = zero IF (arg >= elimit) fn_val = EXP(arg) ELSE ! Use a continued fraction expansion arg = p * LOG(x) - x - lngamma(p) a = one - p b = a + x + one c = zero pn1 = one pn2 = x pn3 = x + one pn4 = x * b fn_val = pn3 / pn4 DO a = a + one b = b + two c = c + one an = a * c pn5 = b * pn3 - an * pn1 pn6 = b * pn4 - an * pn2 IF (ABS(pn6) > zero) THEN rn = pn5 / pn6 IF (ABS(fn_val - rn) <= MIN(tol, tol * rn)) EXIT fn_val = rn END IF pn1 = pn3 pn2 = pn4 pn3 = pn5 pn4 = pn6 IF (ABS(pn5) >= oflo) THEN ! Re-scale terms in continued fraction if terms are large pn1 = pn1 / oflo pn2 = pn2 / oflo pn3 = pn3 / oflo pn4 = pn4 / oflo END IF END DO arg = arg + LOG(fn_val) fn_val = one IF (arg >= elimit) fn_val = one - EXP(arg) END IF RETURN END FUNCTION gammad ! Algorithm AS66 Applied Statistics (1973) vol.22, no.3 ! Evaluates the tail area of the standardised normal curve ! from x to infinity if upper is .true. or ! from minus infinity to x if upper is .false. ! ELF90-compatible version by Alan Miller ! Latest revision - 29 November 2001 FUNCTION alnorm( x, upper ) RESULT( fn_val ) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 100) REAL(dp), INTENT(IN) :: x LOGICAL, INTENT(IN) :: upper REAL(dp) :: fn_val ! Local variables REAL(dp), PARAMETER :: zero=0.0_dp, one=1.0_dp, half=0.5_dp, con=1.28_dp REAL(dp) :: z, y LOGICAL :: up ! Machine dependent constants REAL(dp), PARAMETER :: ltone = 7.0_dp, utzero = 18.66_dp REAL(dp), PARAMETER :: p = 0.398942280444_dp, q = 0.39990348504_dp, & r = 0.398942280385_dp, a1 = 5.75885480458_dp, & a2 = 2.62433121679_dp, a3 = 5.92885724438_dp, & b1 = -29.8213557807_dp, b2 = 48.6959930692_dp, & c1 = -3.8052D-8, c2 = 3.98064794D-4, & c3 = -0.151679116635_dp, c4 = 4.8385912808_dp, & c5 = 0.742380924027_dp, c6 = 3.99019417011_dp, & d1 = 1.00000615302_dp, d2 = 1.98615381364_dp, & d3 = 5.29330324926_dp, d4 = -15.1508972451_dp, & d5 = 30.789933034_dp up = upper z = x IF( z < zero ) THEN up = .NOT. up z = -z END IF IF( z <= ltone .OR. (up .AND. z <= utzero) ) THEN y = half*z*z IF( z > con ) THEN fn_val = r*EXP( -y )/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6)))))) ELSE fn_val = half - z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3)))) END IF ELSE fn_val = zero END IF IF( .NOT. up ) fn_val = one - fn_val RETURN END FUNCTION alnorm FUNCTION lngamma(z) RESULT(lanczos) ! Uses Lanczos-type approximation to ln(gamma) for z > 0. ! Reference: ! Lanczos, C. 'A precision approximation of the gamma ! function', J. SIAM Numer. Anal., B, 1, 86-96, 1964. ! Accuracy: About 14 significant digits except for small regions ! in the vicinity of 1 and 2. ! Programmer: Alan Miller ! 1 Creswick Street, Brighton, Vic. 3187, Australia ! Latest revision - 14 October 1996 IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 60) REAL(dp), INTENT(IN) :: z REAL(dp) :: lanczos ! Local variables REAL(dp) :: a(9) = (/ 0.9999999999995183D0, 676.5203681218835D0, & -1259.139216722289D0, 771.3234287757674D0, & -176.6150291498386D0, 12.50734324009056D0, & -0.1385710331296526D0, 0.9934937113930748D-05, & 0.1659470187408462D-06 /), zero = 0.D0, & one = 1.d0, lnsqrt2pi = 0.9189385332046727D0, & half = 0.5d0, sixpt5 = 6.5d0, seven = 7.d0, tmp INTEGER :: j IF (z <= zero) THEN WRITE(*, *) 'Error: zero or -ve argument for lngamma' RETURN END IF lanczos = zero tmp = z + seven DO j = 9, 2, -1 lanczos = lanczos + a(j)/tmp tmp = tmp - one END DO lanczos = lanczos + a(1) lanczos = LOG(lanczos) + lnsqrt2pi - (z + sixpt5) + (z - half)*LOG(z + sixpt5) RETURN END FUNCTION lngamma