#include "xlisp.h" #include "xlstat.h" #include "linalg.h" #define PLIMIT 1e4 #define XBIG 1e8 #define TOL 1e-14 #define OFLO 1e37 /* * ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3 * Computation of the Incomplete Gamma Integral * Translated by f2c and modified. */ VOID gammabase P3C(double *, x, double *, p, double *, val) { /* Local variables */ double a, b, c, an, rn; double pn1, pn2, pn3, pn4, pn5, pn6, arg; *val = 0.0; if (*x <= 0 || *p <= 0.0) { return; } else if (*p > PLIMIT) { /* Use a normal approximation if P > PLIMIT */ pn1 = sqrt(*p) * 3.0 * (pow(*x / *p, 1.0 / 3.0) + 1.0 / (*p * 9.) - 1.0); normbase(&pn1, val); return; } else if (*x > XBIG) { /* If X is extremely large compared to P then set value = 1 */ *val = 1.0; } else if (*x <= 1.0 || *x < *p) { /* Use Pearson's series expansion. */ /* (Note that P is not large enough to force overflow in gamma). */ arg = *p * log(*x) - *x - gamma(*p + 1.0); c = 1.0; *val = 1.0; a = *p; do { a += 1.0; c = c * *x / a; *val += c; } while (c > TOL); arg += log(*val); *val = exp(arg); } else { /* Use a continued fraction expansion */ arg = *p * log(*x) - *x - gamma(*p); a = 1.0 - *p; b = a + *x + 1.0; c = 0.0; pn1 = 1.0; pn2 = *x; pn3 = *x + 1.0; pn4 = *x * b; *val = pn3 / pn4; while (TRUE) { a += 1.0; b += 2.0; c += 1.0; an = a * c; pn5 = b * pn3 - an * pn1; pn6 = b * pn4 - an * pn2; if (abs(pn6) > 0.0) { rn = pn5 / pn6; if (abs(*val - rn) <= TOL * min(1.0,rn)) break; *val = rn; } pn1 = pn3; pn2 = pn4; pn3 = pn5; pn4 = pn6; if (abs(pn5) >= OFLO) { /* Re-scale terms in continued fraction if terms are large */ pn1 /= OFLO; pn2 /= OFLO; pn3 /= OFLO; pn4 /= OFLO; } } arg += log(*val); *val = 1.0 - exp(arg); } } static double gammad_ P3C(double *, x, double *, a, int *,iflag) { double cdf; gammabase(x, a, &cdf); return(cdf); } /* ppchi2.f -- translated by f2c and modified 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, (but I am using it for 0 < p < 1 - seems to work) v must be positive, g must be supplied and should be equal to ln(gamma(v/2.0)) Auxiliary routines required: ppnd = AS 111 (or AS 241) and gammad_. */ static double ppchi2 P4C(double *, p, double *, v, double *, g, int *, ifault) { /* Initialized data */ static double aa = .6931471806; static double six = 6.; static double c1 = .01; static double c2 = .222222; static double c3 = .32; static double c4 = .4; static double c5 = 1.24; static double c6 = 2.2; static double c7 = 4.67; static double c8 = 6.66; static double c9 = 6.73; static double e = 5e-7; static double c10 = 13.32; static double c11 = 60.; static double c12 = 70.; static double c13 = 84.; static double c14 = 105.; static double c15 = 120.; static double c16 = 127.; static double c17 = 140.; static double c18 = 1175.; static double c19 = 210.; static double c20 = 252.; static double c21 = 2264.; static double c22 = 294.; static double c23 = 346.; static double c24 = 420.; static double c25 = 462.; static double c26 = 606.; static double c27 = 672.; static double c28 = 707.; static double c29 = 735.; static double c30 = 889.; static double c31 = 932.; static double c32 = 966.; static double c33 = 1141.; static double c34 = 1182.; static double c35 = 1278.; static double c36 = 1740.; static double c37 = 2520.; static double c38 = 5040.; static double zero = 0.; static double half = .5; static double one = 1.; static double two = 2.; static double three = 3.; /* static double pmin = 2e-6; static double pmax = .999998; */ static double pmin = 0.0; static double pmax = 1.0; /* System generated locals */ double ret_val, d_1, d_2; /* Local variables */ static double a, b, c, q, t, x, p1, p2, s1, s2, s3, s4, s5, s6, ch; static double xx; static int if1; /* test arguments and initialise */ ret_val = -one; *ifault = 1; if (*p <= pmin || *p >= pmax) return ret_val; *ifault = 2; if (*v <= zero) return ret_val; *ifault = 0; xx = half * *v; c = xx - one; if (*v < -c5 * log(*p)) { /* starting approximation for small chi-squared */ ch = pow(*p * xx * exp(*g + xx * aa), one / xx); if (ch < e) { ret_val = ch; return ret_val; } } else if (*v > c3) { /* call to algorithm AS 111 - note that p has been tested above. */ /* AS 241 could be used as an alternative. */ x = ppnd(*p, &if1); /* starting approximation using Wilson and Hilferty estimate */ p1 = c2 / *v; /* Computing 3rd power */ d_1 = x * sqrt(p1) + one - p1; ch = *v * (d_1 * d_1 * d_1); /* starting approximation for p tending to 1 */ if (ch > c6 * *v + six) ch = -two * (log(one - *p) - c * log(half * ch) + *g); } else{ /* starting approximation for v less than or equal to 0.32 */ ch = c4; a = log(one - *p); do { q = ch; p1 = one + ch * (c7 + ch); p2 = ch * (c9 + ch * (c8 + ch)); d_1 = -half + (c7 + two * ch) / p1; d_2 = (c9 + ch * (c10 + three * ch)) / p2; t = d_1 - d_2; ch -= (one - exp(a + *g + half * ch + c * aa) * p2 / p1) / t; } while (fabs(q / ch - one) > c1); } do { /* call to gammad_ and calculation of seven term Taylor series */ q = ch; p1 = half * ch; p2 = *p - gammad_(&p1, &xx, &if1); if (if1 != 0) { *ifault = 3; return ret_val; } t = p2 * exp(xx * aa + *g + p1 - c * log(ch)); b = t / ch; 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; d_1 = (s3 - b * (s4 - b * (s5 - b * s6))); d_1 = (s1 - b * (s2 - b * d_1)); ch += t * (one + half * t * s1 - b * c * d_1); } while (fabs(q / ch - one) > e); ret_val = ch; return ret_val; } double ppgamma P3C(double, p, double, a, int *, ifault) { double x, v, g; if (p == 0.0) return 0.0; else { v = 2.0 * a; g = gamma(a); x = ppchi2(&p, &v, &g, ifault); return x / 2.0; } }