lib/gammabase.c

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