src/c/gammab.c

   1 #include "xlisp.h"
   2 #include "xlstat.h"
   3 #include "linalg.h"
   4
   5 #define PLIMIT 1e4
   6 #define XBIG 1e8
   7 #define TOL 1e-14
   8 #define OFLO 1e37
   9
  10 /*
  11  * ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3
  12  * Computation of the Incomplete Gamma Integral
  13  * Translated by f2c and modified.
  14  */
  15
  16 VOID gammabase P3C(double *, x, double *, p, double *, val)
  17 {
  18   /* Local variables */
  19   double a, b, c, an, rn;
  20   double pn1, pn2, pn3, pn4, pn5, pn6, arg;
  21
  22   *val = 0.0;
  23   if (*x <= 0 || *p <= 0.0) {
  24     return;
  25   }
  26   else if (*p > PLIMIT) {
  27     /* Use a normal approximation if P > PLIMIT */
  28     pn1 = sqrt(*p) * 3.0 * (pow(*x / *p, 1.0 / 3.0) + 1.0 / (*p * 9.) - 1.0);
  29     normbase(&pn1, val);
  30     return;
  31   }
  32   else if (*x > XBIG) {
  33     /* If X is extremely large compared to P then set value = 1 */
  34     *val = 1.0;
  35   }
  36   else if (*x <= 1.0 || *x < *p) {
  37     /* Use Pearson's series expansion. */
  38     /* (Note that P is not large enough to force overflow in gamma). */
  39     arg = *p * log(*x) - *x - gamma(*p + 1.0);
  40     c = 1.0;
  41     *val = 1.0;
  42     a = *p;
  43     do {
  44       a += 1.0;
  45       c = c * *x / a;
  46       *val += c;
  47     } while (c > TOL);
  48     arg += log(*val);
  49     *val = exp(arg);
  50   }
  51   else {
  52     /* Use a continued fraction expansion */
  53     arg = *p * log(*x) - *x - gamma(*p);
  54     a = 1.0 - *p;
  55     b = a + *x + 1.0;
  56     c = 0.0;
  57     pn1 = 1.0;
  58     pn2 = *x;
  59     pn3 = *x + 1.0;
  60     pn4 = *x * b;
  61     *val = pn3 / pn4;
  62     while (TRUE) {
  63       a += 1.0;
  64       b += 2.0;
  65       c += 1.0;
  66       an = a * c;
  67       pn5 = b * pn3 - an * pn1;
  68       pn6 = b * pn4 - an * pn2;
  69       if (abs(pn6) > 0.0) {
  70         rn = pn5 / pn6;
  71         if (abs(*val - rn) <= TOL * min(1.0,rn))
  72           break;
  73         *val = rn;
  74       }
  75
  76       pn1 = pn3;
  77       pn2 = pn4;
  78       pn3 = pn5;
  79       pn4 = pn6;
  80       if (abs(pn5) >= OFLO) {
  81         /* Re-scale terms in continued fraction if terms are large */
  82         pn1 /= OFLO;
  83         pn2 /= OFLO;
  84         pn3 /= OFLO;
  85         pn4 /= OFLO;
  86       }
  87     }
  88     arg += log(*val);
  89     *val = 1.0 - exp(arg);
  90   }
  91 }
  92
  93 static double gammad_ P3C(double *, x, double *, a, int *,iflag)
  94 {
  95   double cdf;
  96
  97   gammabase(x, a, &cdf);
  98   return(cdf);
  99 }
 100
 101 /*
 102   ppchi2.f -- translated by f2c and modified
 103
 104   Algorithm AS 91   Appl. Statist. (1975) Vol.24, P.35
 105   To evaluate the percentage points of the chi-squared
 106   probability distribution function.
 107
 108   p must lie in the range 0.000002 to 0.999998,
 109     (but I am using it for 0 < p < 1 - seems to work)
 110   v must be positive,
 111   g must be supplied and should be equal to ln(gamma(v/2.0))
 112
 113   Auxiliary routines required: ppnd = AS 111 (or AS 241) and gammad_.
 114 */
 115
 116 static double ppchi2 P4C(double *, p, double *, v, double *, g, int *, ifault)
 117 {
 118   /* Initialized data */
 119
 120   static double aa = .6931471806;
 121   static double six = 6.;
 122   static double c1 = .01;
 123   static double c2 = .222222;
 124   static double c3 = .32;
 125   static double c4 = .4;
 126   static double c5 = 1.24;
 127   static double c6 = 2.2;
 128   static double c7 = 4.67;
 129   static double c8 = 6.66;
 130   static double c9 = 6.73;
 131   static double e = 5e-7;
 132   static double c10 = 13.32;
 133   static double c11 = 60.;
 134   static double c12 = 70.;
 135   static double c13 = 84.;
 136   static double c14 = 105.;
 137   static double c15 = 120.;
 138   static double c16 = 127.;
 139   static double c17 = 140.;
 140   static double c18 = 1175.;
 141   static double c19 = 210.;
 142   static double c20 = 252.;
 143   static double c21 = 2264.;
 144   static double c22 = 294.;
 145   static double c23 = 346.;
 146   static double c24 = 420.;
 147   static double c25 = 462.;
 148   static double c26 = 606.;
 149   static double c27 = 672.;
 150   static double c28 = 707.;
 151   static double c29 = 735.;
 152   static double c30 = 889.;
 153   static double c31 = 932.;
 154   static double c32 = 966.;
 155   static double c33 = 1141.;
 156   static double c34 = 1182.;
 157   static double c35 = 1278.;
 158   static double c36 = 1740.;
 159   static double c37 = 2520.;
 160   static double c38 = 5040.;
 161   static double zero = 0.;
 162   static double half = .5;
 163   static double one = 1.;
 164   static double two = 2.;
 165   static double three = 3.;
 166
 167 /*
 168   static double pmin = 2e-6;
 169   static double pmax = .999998;
 170 */
 171   static double pmin = 0.0;
 172   static double pmax = 1.0;
 173
 174   /* System generated locals */
 175   double ret_val, d_1, d_2;
 176
 177   /* Local variables */
 178   static double a, b, c, q, t, x, p1, p2, s1, s2, s3, s4, s5, s6, ch;
 179   static double xx;
 180   static int if1;
 181
 182
 183   /* test arguments and initialise */
 184   ret_val = -one;
 185   *ifault = 1;
 186   if (*p <= pmin || *p >= pmax) return ret_val;
 187   *ifault = 2;
 188   if (*v <= zero) return ret_val;
 189   *ifault = 0;
 190   xx = half * *v;
 191   c = xx - one;
 192
 193   if (*v < -c5 * log(*p)) {
 194     /* starting approximation for small chi-squared */
 195     ch = pow(*p * xx * exp(*g + xx * aa), one / xx);
 196     if (ch < e) {
 197       ret_val = ch;
 198       return ret_val;
 199     }
 200   }
 201   else if (*v > c3) {
 202     /* call to algorithm AS 111 - note that p has been tested above. */
 203     /* AS 241 could be used as an alternative. */
 204     x = ppnd(*p, &if1);
 205
 206     /* starting approximation using Wilson and Hilferty estimate */
 207     p1 = c2 / *v;
 208     /* Computing 3rd power */
 209     d_1 = x * sqrt(p1) + one - p1;
 210     ch = *v * (d_1 * d_1 * d_1);
 211
 212     /* starting approximation for p tending to 1 */
 213     if (ch > c6 * *v + six)
 214       ch = -two * (log(one - *p) - c * log(half * ch) + *g);
 215   }
 216   else{
 217     /* starting approximation for v less than or equal to 0.32 */
 218     ch = c4;
 219     a = log(one - *p);
 220     do {
 221       q = ch;
 222       p1 = one + ch * (c7 + ch);
 223       p2 = ch * (c9 + ch * (c8 + ch));
 224       d_1 = -half + (c7 + two * ch) / p1;
 225       d_2 = (c9 + ch * (c10 + three * ch)) / p2;
 226       t = d_1 - d_2;
 227       ch -= (one - exp(a + *g + half * ch + c * aa) * p2 / p1) / t;
 228     } while (fabs(q / ch - one) > c1);
 229   }
 230
 231   do {
 232     /* call to gammad_ and calculation of seven term Taylor series */
 233     q = ch;
 234     p1 = half * ch;
 235     p2 = *p - gammad_(&p1, &xx, &if1);
 236     if (if1 != 0) {
 237       *ifault = 3;
 238       return ret_val;
 239     }
 240     t = p2 * exp(xx * aa + *g + p1 - c * log(ch));
 241     b = t / ch;
 242     a = half * t - b * c;
 243     s1 = (c19 + a * (c17 + a * (c14 + a * (c13 + a * (c12 + c11 * a))))) / c24;
 244     s2 = (c24 + a * (c29 + a * (c32 + a * (c33 + c35 * a)))) / c37;
 245     s3 = (c19 + a * (c25 + a * (c28 + c31 * a))) / c37;
 246     s4 = (c20 + a * (c27 + c34 * a) + c * (c22 + a * (c30 + c36 * a))) / c38;
 247     s5 = (c13 + c21 * a + c * (c18 + c26 * a)) / c37;
 248     s6 = (c15 + c * (c23 + c16 * c)) / c38;
 249     d_1 = (s3 - b * (s4 - b * (s5 - b * s6)));
 250     d_1 = (s1 - b * (s2 - b * d_1));
 251     ch += t * (one + half * t * s1 - b * c * d_1);
 252   } while (fabs(q / ch - one) > e);
 253
 254   ret_val = ch;
 255   return ret_val;
 256 }
 257
 258 double ppgamma P3C(double, p, double, a, int *, ifault)
 259 {
 260   double x, v, g;
 261
 262   if (p == 0.0)
 263     return 0.0;
 264   else {
 265     v = 2.0 * a;
 266     g = gamma(a);
 267     x = ppchi2(&p, &v, &g, ifault);
 268     return x / 2.0;
 269   }
 270 }