/*************************************************************************/ /* */ /* ALGORITHM AS 63: THE INCOMPLETE BETA INTEGRAL */ /* */ /* K. L. Majunder & G. P. Bhattacharjee */ /* Department of Mathematics, I.I.T., Kharagpur, India. */ /* */ /*************************************************************************/ /* Journal of the Royal Statistical Society (Series C): */ /* Applied Statistics, Vol.22, No. 3, pp. 409-411 (1973). */ /* Log of the beta function */ /* Includes log of the gamma function */ betaln: PROCEDURE (a0, b0) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* ----------------------------------------------------------------------- */ /* EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION */ /* ----------------------------------------------------------------------- */ /* E = 0.5*LN(2*PI) */ /* -------------------------- */ DECLARE ( a0, b0 ) FLOAT (18); /* Local variables */ DECLARE ( e STATIC INITIAL ( .918938533204673E0) ) FLOAT (18); DECLARE ( a, b, c, h, u, v, w, z ) FLOAT (18); DECLARE ( i, n ) FIXED BINARY (31); DECLARE ( half STATIC INITIAL ( 0.5), one STATIC INITIAL ( 1.0), two STATIC INITIAL ( 2.0), eight STATIC INITIAL ( 8.0) ) FLOAT (18); /* -------------------------- */ a = MIN(a0, b0); b = MAX(a0, b0); IF a < eight THEN DO; IF a < one THEN DO; /* ----------------------------------------------------------------------- */ /* PROCEDURE WHEN A < 1 */ /* ----------------------------------------------------------------------- */ IF b < eight THEN RETURN ( gamln(a) + (gamln(b) - gamln(a+b)) ); RETURN ( gamln(a) + algdiv(a,b) ); END; /* ----------------------------------------------------------------------- */ /* PROCEDURE WHEN 1 <= A < 8 */ /* ----------------------------------------------------------------------- */ IF a <= two THEN DO; IF b <= two THEN RETURN (gamln(a) + gamln(b) - gsumln(a,b) ); w = 0; IF b < eight THEN GO TO L20; RETURN (gamln(a) + algdiv(a,b)); END; /* REDUCTION OF A WHEN B <= 1000 */ IF b > 1000.00000000000000E0 THEN GO TO L40; n = a - one; w = one; DO i = 1 TO n; a = a - one; h = a / b; w = w * (h/(one + h)); END; w = LOG(w); IF b >= eight THEN RETURN ( w + gamln(a) + algdiv(a,b) ); /* REDUCTION OF B WHEN B < 8 */ L20: n = b - one; z = one; DO i = 1 TO n; b = b - one; z = z * (b/(a+b)); END; RETURN ( w + LOG(z) + (gamln(a) + (gamln(b) - gsumln(a,b))) ); /* REDUCTION OF A WHEN B > 1000 */ L40: n = a - one; w = one; DO i = 1 TO n; a = a - one; w = w * (a/(one+a/b)); END; RETURN ( (LOG(w) - n*LOG(b)) + (gamln(a) + algdiv(a,b)) ); END; /* ----------------------------------------------------------------------- */ /* PROCEDURE WHEN A >= 8 */ /* ----------------------------------------------------------------------- */ w = bcorr(a,b); h = a / b; c = h / (one + h); u = -(a-half) * LOG(c); v = b * alnrel(h); IF u > v THEN RETURN ( (((-half*LOG(b) + e) + w) - v) - u ); RETURN ( (((-half*LOG(b) + e) + w) - u) - v ); END betaln; gamln: PROCEDURE (a) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* ----------------------------------------------------------------------- */ /* EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A */ /* ----------------------------------------------------------------------- */ /* WRITTEN BY ALFRED H. MORRIS */ /* NAVAL SURFACE WARFARE CENTER */ /* DAHLGREN, VIRGINIA */ /* -------------------------- */ /* D = 0.5*(LN(2*PI) - 1) */ /* -------------------------- */ DECLARE ( a ) FLOAT (18); DECLARE ( fn_val ) FLOAT (18); /* Local variables */ DECLARE ( d STATIC INITIAL ( 0.418938533204673E0) , c0 STATIC INITIAL ( 0.833333333333333E-01) , c1 STATIC INITIAL ( -0.277777777760991E-02) , c2 STATIC INITIAL ( 0.793650666825390E-03) , c3 STATIC INITIAL ( -0.595202931351870E-03) , c4 STATIC INITIAL ( 0.837308034031215E-03) , c5 STATIC INITIAL ( -0.165322962780713E-02) ) FLOAT (18); DECLARE ( t, w ) FLOAT (18); DECLARE ( i, n ) FIXED BINARY (31); /* ----------------------------------------------------------------------- */ IF a <= 0.80000000000000000E0 THEN RETURN ( gamln1(a) - LOG(a) ); IF a <= 2.25000000000000000E0 THEN DO; t = (a-0.50000000000000000E0) - 0.50000000000000000E0; RETURN ( gamln1(t) ); END; IF a < 10.0000000000000000E0 THEN DO; n = a - 1.25000000000000000E0; t = a; w = 1.00000000000000000E0; DO i = 1 TO n; t = t - 1.00000000000000000E0; w = t * w; END; RETURN ( gamln1(t-1.00000000000000000E0) + LOG(w) ); END; t = (1.00000000000000000E0/a) ** 2; w = (((((c5*t + c4)*t + c3)*t + c2)*t + c1)*t + c0) / a; RETURN ( (d+w) + (a-0.50000000000000000E0) * (LOG(a)-1.00000000000000000E0) ); END gamln; algdiv: PROCEDURE (a, b) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* ----------------------------------------------------------------------- */ /* COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B >= 8 */ /* -------- */ /* IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY */ /* LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X). */ /* ----------------------------------------------------------------------- */ DECLARE ( a, b ) FLOAT (18); DECLARE ( fn_val ) FLOAT (18); /* Local variables */ DECLARE ( c0 STATIC INITIAL ( .833333333333333E-01) , c1 STATIC INITIAL ( -.277777777760991E-02) , c2 STATIC INITIAL ( .793650666825390E-03) , c3 STATIC INITIAL ( -.595202931351870E-03) , c4 STATIC INITIAL ( .837308034031215E-03) , c5 STATIC INITIAL ( -.165322962780713E-02) ) FLOAT (18); DECLARE ( c, d, h, s3, s5, s7, s9, s11, t, u, v, w, x, x2 ) FLOAT (18); /* ------------------------ */ IF a > b THEN DO; h = b / a; c = 1.00000000000000000E0 / (1.00000000000000000E0 + h); x = h / (1.00000000000000000E0 + h); d = a + (b - 0.50000000000000000E0); END; ELSE DO; h = a / b; c = h / (1.00000000000000000E0 + h); x = 1.00000000000000000E0 / (1.00000000000000000E0 + h); d = b + (a - 0.50000000000000000E0); END; /* SET SN = (1 - X**N)/(1 - X) */ x2 = x * x; s3 = 1.00000000000000000E0 + (x + x2); s5 = 1.00000000000000000E0 + (x + x2*s3); s7 = 1.00000000000000000E0 + (x + x2*s5); s9 = 1.00000000000000000E0 + (x + x2*s7); s11 = 1.00000000000000000E0 + (x + x2*s9); /* SET W = DEL(B) - DEL(A + B) */ t = (1.00000000000000000E0/b) ** 2; w = ((((c5*s11*t + c4*s9)*t + c3*s7)*t + c2*s5)*t + c1*s3) * t + c0; w = w * (c/b); /* COMBINE THE RESULTS */ u = d * alnrel(a/b); v = a * (LOG(b) - 1.00000000000000000E0); IF u > v THEN RETURN ( (w-v) - u ); RETURN ( (w-u) - v ); END algdiv; gsumln: PROCEDURE (a, b) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* ----------------------------------------------------------------------- */ /* EVALUATION OF THE FUNCTION LN(GAMMA(A + B)) */ /* FOR 1 <= A <= 2 AND 1 <= B <= 2 */ /* ----------------------------------------------------------------------- */ DECLARE ( a, b ) FLOAT (18); DECLARE ( fn_val ) FLOAT (18); /* Local variables */ DECLARE ( x ) FLOAT (18); x = a + b - 2.00000000000000000E0; IF x <= 0.25000000000000000E0 THEN RETURN ( gamln1(1.00000000000000000E0 + x) ); IF x <= 1.25000000000000000E0 THEN RETURN ( gamln1(x) + alnrel(x) ); RETURN ( gamln1(x - 1.00000000000000000E0) + LOG(x*(1.00000000000000000E0 + x)) ); END gsumln; bcorr: PROCEDURE (a0, b0) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* ----------------------------------------------------------------------- */ /* EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE */ /* LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A). */ /* IT IS ASSUMED THAT A0 >= 8 AND B0 >= 8. */ /* ----------------------------------------------------------------------- */ DECLARE ( a0, b0 ) FLOAT (18); /* Local variables */ DECLARE ( c0 STATIC INITIAL ( .833333333333333E-01) , c1 STATIC INITIAL ( -.277777777760991E-02) , c2 STATIC INITIAL ( .793650666825390E-03) , c3 STATIC INITIAL ( -.595202931351870E-03) , c4 STATIC INITIAL ( .837308034031215E-03) , c5 STATIC INITIAL ( -.165322962780713E-02) ) FLOAT (18); DECLARE ( a, b, c, h, s3, s5, s7, s9, s11, t, w, x, x2 ) FLOAT (18); /* ------------------------ */ a = MIN(a0,b0); b = MAX(a0,b0); h = a / b; c = h / (1.00000000000000000E0 + h); x = 1.00000000000000000E0 / (1.00000000000000000E0 + h); x2 = x * x; /* SET SN = (1 - X**N)/(1 - X) */ s3 = 1.00000000000000000E0 + (x + x2); s5 = 1.00000000000000000E0 + (x + x2*s3); s7 = 1.00000000000000000E0 + (x + x2*s5); s9 = 1.00000000000000000E0 + (x + x2*s7); s11 = 1.00000000000000000E0 + (x + x2*s9); /* SET W = DEL(B) - DEL(A + B) */ t = (1.00000000000000000E0/b) ** 2; w = ((((c5*s11*t + c4*s9)*t + c3*s7)*t + c2*s5)*t + c1*s3) * t + c0; w = w * (c/b); /* COMPUTE DEL(A) + W */ t = (1.00000000000000000E0/a) ** 2; RETURN ( (((((c5*t + c4)*t + c3)*t + c2)*t + c1)*t + c0) / a + w ); END bcorr; alnrel: PROCEDURE (a) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* ----------------------------------------------------------------------- */ /* EVALUATION OF THE FUNCTION LN(1 + A) */ /* ----------------------------------------------------------------------- */ DECLARE ( a ) FLOAT (18); /* Local variables */ DECLARE ( p1 STATIC INITIAL ( -.129418923021993E+01) , p2 STATIC INITIAL ( .405303492862024E+00) , p3 STATIC INITIAL ( -.178874546012214E-01) , q1 STATIC INITIAL ( -.162752256355323E+01) , q2 STATIC INITIAL ( .747811014037616E+00) , q3 STATIC INITIAL ( -.845104217945565E-01) , zero STATIC INITIAL ( 0.0) , half STATIC INITIAL ( 0.5) , one STATIC INITIAL ( 1.0) , two STATIC INITIAL ( 2.0) ) FLOAT (18); DECLARE ( t, t2, w, x ) FLOAT (18); /* -------------------------- */ IF ABS(a) <= 0.37500000000000000E0 THEN DO; t = a / (a + two); t2 = t * t; w = (((p3*t2 + p2)*t2 + p1)*t2 + one) / (((q3*t2 + q2)*t2 + q1)*t2 + one); RETURN ( two * t * w ); END; x = one + a; IF a < zero THEN x = (a + half) + half; RETURN ( LOG(x) ); END alnrel; gamln1: PROCEDURE (a) RETURNS( FLOAT (18) ) OPTIONS (REORDER); /* ----------------------------------------------------------------------- */ /* EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 <= A <= 1.25 */ /* ----------------------------------------------------------------------- */ DECLARE ( a ) FLOAT (18); /* Local variables */ DECLARE ( p0 STATIC INITIAL ( .577215664901533E+00) , p1 STATIC INITIAL ( .844203922187225E+00) , p2 STATIC INITIAL ( -.168860593646662E+00) , p3 STATIC INITIAL ( -.780427615533591E+00) , p4 STATIC INITIAL ( -.402055799310489E+00) , p5 STATIC INITIAL ( -.673562214325671E-01) , p6 STATIC INITIAL ( -.271935708322958E-02) , q1 STATIC INITIAL ( .288743195473681E+01) , q2 STATIC INITIAL ( .312755088914843E+01) , q3 STATIC INITIAL ( .156875193295039E+01) , q4 STATIC INITIAL ( .361951990101499E+00) , q5 STATIC INITIAL ( .325038868253937E-01) , q6 STATIC INITIAL ( .667465618796164E-03) , r0 STATIC INITIAL ( .422784335098467E+00) , r1 STATIC INITIAL ( .848044614534529E+00) , r2 STATIC INITIAL ( .565221050691933E+00) , r3 STATIC INITIAL ( .156513060486551E+00) , r4 STATIC INITIAL ( .170502484022650E-01) , r5 STATIC INITIAL ( .497958207639485E-03) , s1 STATIC INITIAL ( .124313399877507E+01) , s2 STATIC INITIAL ( .548042109832463E+00) , s3 STATIC INITIAL ( .101552187439830E+00) , s4 STATIC INITIAL ( .713309612391000E-02) , s5 STATIC INITIAL ( .116165475989616E-03) ) FLOAT (18); DECLARE ( w, x ) FLOAT (18); /* ---------------------- */ IF a < 0.60000000000000000E0 THEN DO; w = ((((((p6*a + p5)*a + p4)*a + p3)*a + p2)*a + p1)*a + p0) / ((((((q6*a + q5)*a + q4)*a + q3)*a + q2)*a + q1)*a + 1.00000000000000000E0); RETURN ( -a * w ); END; x = (a - 0.50000000000000000E0) - 0.50000000000000000E0; w = (((((r5*x + r4)*x + r3)*x + r2)*x + r1)*x + r0) / (((((s5*x + s4)*x + s3)*x + s2)*x + s1)*x + 1.00000000000000000E0); RETURN ( x * w ); END gamln1; betain: PROCEDURE (x, p, q, beta) RETURNS( FLOAT (18) ); /* Algorithm AS 63 Appl. Statist. (1973), vol.22, no.3 */ /* Computes incomplete beta function ratio for arguments */ /* x between zero and one, p and q positive. */ /* Log of complete beta function, beta, is assumed to be known */ /* N.B. Argument IFAULT has been removed */ /* Latest revision - 5 July 2003 */ DECLARE ( x, p, q, beta ) FLOAT (18); DECLARE ( fn_val ) FLOAT (18); /* Local variables */ DECLARE ( indx ) BIT(1) ALIGNED; DECLARE ( ns ) FIXED BINARY (31); DECLARE ( psq, cx, xx, pp, qq, term, ai, rx, temp ) FLOAT (18); /* Define accuracy and initialise */ DECLARE ( zero STATIC INITIAL ( 0.0) , one STATIC INITIAL ( 1.0) , acu STATIC INITIAL ( 1.0E-14) ) FLOAT (18); fn_val = x; /* Test for admissibility of arguments */ IF p <= zero | q <= zero THEN DO; PUT SKIP LIST ('AS63: Either p or q <= 0'); RETURN (x); END; IF x < zero | x > one THEN DO; PUT SKIP LIST ('AS63: Argument x outside range (0, 1)'); RETURN (x); END; IF x = zero | x = one THEN RETURN (x); /* Change tail if necessary and determine s */ psq = p + q; cx = one - x; IF p < psq*x THEN DO; xx = cx; cx = x; pp = q; qq = p; indx = '1'B; END; ELSE DO; xx = x; pp = p; qq = q; indx = '0'B; END; term = one; ai = one; fn_val = one; ns = qq + cx*psq; /* Use Soper's reduction formulae. */ rx = xx/cx; L3: temp = qq - ai; IF ns = 0 THEN rx = xx; L4: term = term*temp*rx/(pp+ai); fn_val = fn_val + term; temp = ABS(term); IF temp <= acu & temp <= acu*fn_val THEN GO TO L5; ai = ai + one; ns = ns - 1; IF ns >= 0 THEN GO TO L3; temp = psq; psq = psq+one; GO TO L4; /* Calculate result */ L5: fn_val = fn_val*EXP(pp*LOG(xx) + (qq-one)*LOG(cx) - beta)/pp; IF indx THEN RETURN (one - fn_val); RETURN (fn_val); END betain;