/*************************************************************************/ /* */ /* ALGORITHM AS239 CHI-SQUARED AND INCOMPLETE GAMMA INTEGRAL */ /* */ /* B. L. Shea, Manchester Polytechnis, U.K. */ /* */ /*************************************************************************/ /* Journal of the Royal Statistical Society (Series C): */ /* Applied Statistics, Vol.37 No. 3, 1988, pp. 466-472. */ /* Computation of the Incomplete Gamma Integral */ /* Translated from Fortran 90 to PL/I by R. A. Vowels, 3 April 2006. */ /* ELF90-compatible version by Alan Miller */ /* Latest revision - 27 October 2000 */ gammad: PROCEDURE (x, p) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* N.B. Argument IFAULT has been removed */ /* It was, but it is still used, and wasn't made a local variable */ /* (Translator's note, 23 May 2007). */ /* Auxiliary functions required: ALOGAM = logarithm of the gamma */ /* function, and ALNORM = algorithm AS66 */ DECLARE ( x, p ) FLOAT (18); DECLARE ( fn_val ) FLOAT (18); /* Local variables */ DECLARE ifault FIXED BINARY; DECLARE ( pn1, pn2, pn3, pn4, pn5, pn6, arg, c, rn, a, b, an ) FLOAT (18); DECLARE ( zero STATIC INITIAL ( 0.00000000000000000E0) , one STATIC INITIAL ( 1.00000000000000000E0) , two STATIC INITIAL ( 2.00000000000000000E0) , oflo STATIC INITIAL ( 1.00000000000000000E+37) , three STATIC INITIAL ( 3.00000000000000000E0) , nine STATIC INITIAL ( 9.00000000000000000E0) , tol STATIC INITIAL ( 1.00000000000000000E-14) , xbig STATIC INITIAL ( 1.00000000000000000E+8) , plimit STATIC INITIAL ( 1000.00000000000000E0) , elimit STATIC INITIAL ( -88.0000000000000000E0) ) FLOAT (18); /* EXTERNAL alogam, alnorm */ fn_val = zero; /* Check that we have valid values for X and P */ IF p <= zero | x < zero THEN DO; PUT SKIP LIST ('AS 239: Either p <= 0 or x < 0'); RETURN (fn_val); END; IF x = zero THEN RETURN (fn_val); /* Use a normal approximation if P > PLIMIT */ IF p > plimit THEN DO; pn1 = three * SQRT(p) * ((x / p) ** (one / three) + one /(nine * p) - one); RETURN ( alnorm(pn1, '0'B ) ); END; /* If X is extremely large compared to P then set fn_val = 1 */ IF x > xbig THEN RETURN (one); IF x <= one | x < p THEN DO; /* Use Pearson's series expansion. */ /* (Note that P is not large enough to force overflow in ALOGAM). */ /* No need to test IFAULT on exit since P > 0. */ arg = p * LOG(x) - x - alogam(p + one, ifault); c = one; fn_val = one; a = p; L40: a = a + one; c = c * x / a; fn_val = fn_val + c; IF c > tol THEN GO TO L40; arg = arg + LOG(fn_val); fn_val = zero; IF arg >= elimit THEN fn_val = EXP(arg); END; ELSE DO; /* Use a continued fraction expansion */ arg = p * LOG(x) - x - alogam(p, ifault); a = one - p; b = a + x + one; c = zero; pn1 = one; pn2 = x; pn3 = x + one; pn4 = x * b; fn_val = pn3 / pn4; L60: 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 DO; rn = pn5 / pn6; IF ABS(fn_val - rn) <= MIN(tol, tol * rn) THEN GO TO L80; fn_val = rn; END; pn1 = pn3; pn2 = pn4; pn3 = pn5; pn4 = pn6; IF ABS(pn5) >= oflo THEN DO; /* Re-scale terms in continued fraction if terms are large */ pn1 = pn1 / oflo; pn2 = pn2 / oflo; pn3 = pn3 / oflo; pn4 = pn4 / oflo; END; GO TO L60; L80: arg = arg + LOG(fn_val); fn_val = one; IF arg >= elimit THEN fn_val = one - EXP(arg); END; RETURN (fn_val); END gammad;