lib/studentbase.c

   1 #include "xmath.h"
   2
   3 #define TWOVRPI 0.636619772367581343
   4 #define HALF_PI 1.5707963268
   5 #define EPSILON .000001
   6
   7 extern void betabase(double *, double *, double *,int *, int *,double *);
   8 extern void normbase(double *,double *);
   9
  10 extern double ppnd(), ppbeta();
  11
  12 /* CACM Algorithm 395, by G. W. Hill */
  13
  14 void
  15 studentbase(double *x, double *df, double *cdf)
  16 {
  17   double t, y, b, a, z, j, n;
  18
  19   n = *df;
  20   z = 1.0;
  21   t = (*x) * (*x);
  22   y = t / n;
  23   b = 1.0 + y;
  24
  25   if (n > floor(n) || (n >= 20 && t < n) || (n > 20)) {
  26     if (n < 2.0 && n != 1.0) {
  27       /* beta integral aproximation for small df */
  28       double  da = 0.5, db = 0.5 * n, dx, dp;
  29       int ia = 0, ib = floor(db);
  30
  31       dx = db / (db + da * t);
  32       betabase(&dx, &db, &da, &ib, &ia, &dp);
  33       *cdf = (*x >= 0) ? 1.0 - .5 * dp : .5 * dp;
  34     } else {
  35       /* asymptotic series for large or non-integer df */
  36       if(y > EPSILON) {
  37         y = log(b);
  38       }
  39       a = n - 0.5;
  40       b = 48.0 * a * a;
  41       y = a * y;
  42       y = (((((-0.4 * y - 3.3) * y - 24.0) * y - 85.5)
  43             / (0.8 * y * y + 100.0 + b) + y + 3.0) / b + 1.0) * sqrt(y);
  44
  45       y = -1.0 * y;
  46       normbase(&y, cdf);
  47       if (*x > 0.0) *cdf = 1.0 - *cdf;
  48     }
  49   } else {
  50     /* nested summation of cosine series */
  51     if (n < 20 && t < 4.0) {
  52       a = y = sqrt(y);
  53       if(n == 1.0) a = 0.0;
  54     } else {
  55       a = sqrt(b);
  56       y = a * n;
  57       for(j = 2; fabs(a - z) > EPSILON; j += 2.0) {
  58         z = a;
  59         y = (y * (j - 1)) / (b * j);
  60         a = a + y / (n + j);
  61       }
  62       n += 2.0;
  63       z = y = 0.0;
  64       a = -a;
  65     }
  66     for(n = n - 2.0; n > 1.0; n -= 2.0) a = ((n - 1.0) / (b * n)) * a + y;
  67     a = (fabs(n) < EPSILON) ? a/sqrt(b) : TWOVRPI * (atan(y) + a / b);
  68     *cdf = z - a;
  69     if(*x > 0.0) {
  70       *cdf = 1.0 - 0.5 * *cdf;
  71     } else {
  72       *cdf = 0.5 * *cdf;
  73     }
  74   }
  75 }
  76
  77 /* CACM Algorithm 396, by G. W. Hill */
  78
  79 double ppstudent(pp, n, ifault)
  80      double pp, n;
  81      int *ifault;
  82 {
  83   double sq, p, a, b, c, d, x, y;
  84
  85   /* convert to double upper tailed probability */
  86   p = (pp < 0.5) ? 2.0 * pp : 2.0 * (1.0 - pp);
  87
  88   if (n <= 3.0) {
  89     if (n == 1) sq = tan(HALF_PI * (1.0 - p));
  90     else if (n == 2.0) sq = sqrt(2.0 / (p * (2.0 - p)) - 2.0);
  91     else {
  92       sq = ppbeta(p, 0.5 * n, 0.5, ifault);
  93       if (sq != 0.0) sq = sqrt((n / sq) - n);
  94     }
  95   }
  96   else {
  97     a = 1.0 / (n - 0.5);
  98     b = 48.0 / (a * a);
  99     c = ((20700.0 * a / b - 98.0) * a - 16) * a + 96.36;
 100     d = ((94.5 / (b + c) - 3.0) / b + 1.0) * sqrt(a * HALF_PI) * n;
 101     x = d * p;
 102     y = pow(x, 2.0 / n);
 103     if (y > 0.05 + a) {
 104       /* asymptotic inverse expansion about normal */
 105       x = ppnd(0.5 * p, ifault);
 106       y = x * x;
 107       if (n < 5) c = c + 0.3 * (n - 4.5) * (x + 0.6);
 108       c = (((0.05 * d * x - 5.0) * x - 7.0) * x - 2.0) * x + b + c;
 109       y = (((((0.4 * y + 6.3) * y + 36.0) * y + 94.5) / c - y - 3.0) / b
 110            + 1.0) * x;
 111       y = a * y * y;
 112       y = (y > .002) ? exp(y) - 1.0 : 0.5 * y * y + y;
 113     }
 114     else {
 115       y = ((1.0 / (((n + 6.0) / (n * y) - 0.089 * d - 0.822)
 116                    * (n + 2.0) * 3.0) + 0.5 / (n + 4.0)) * y - 1.0)
 117         * (n + 1.0) / (n + 2.0) + 1.0 / y;
 118     }
 119     sq = sqrt(n * y);
 120   }
 121
 122   /* decode based on tail */
 123   if (pp < 0.5) sq = -sq;
 124   return(sq);
 125 }