lib/gammabase.c

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