lib/betabase.c

   1 #include "xmath.h"
   2
   3 #define TRUE 1
   4 #define FALSE 0
   5
   6 /* external function */
   7 extern double macheps(), gamma();
   8
   9 /* forward declarations */
  10 static double logbeta(), betai(), xinbta();
  11
  12 void
  13 betabase(double *x, double *a, double *b, int *gia, int *gib, double *cdf)
  14 {
  15   *cdf = betai(*x, *a, *b);
  16 }
  17
  18 double ppbeta(double p, double a, double b, int *ifault)
  19 {
  20   double lbeta;
  21
  22   lbeta = gamma(a) + gamma(b) - gamma(a + b);
  23   return(xinbta(&a, &b, &lbeta, &p, ifault));
  24 }
  25
  26 /*
  27   Static routines
  28 */
  29
  30 static double logbeta(p, q)
  31      double p, q;
  32 {
  33   return(gamma(p) + gamma(q) - gamma(p + q));
  34 }
  35
  36 #define Min(x,y) (((x) < (y)) ? (x) : (y))
  37 #define Max(x,y) (((x) > (y)) ? (x) : (y))
  38
  39 static double betai(x, pin, qin)
  40      double x, pin, qin;
  41 {
  42   /* Translated from FORTRAN
  43      july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
  44      based on bosten and battiste, remark on algorithm 179, comm. acm,
  45      v 17, p 153, (1974).
  46
  47      input arguments --
  48      x      upper limit of integration.  x must be in (0,1) inclusive.
  49      p      first beta distribution parameter.  p must be gt 0.0.
  50      q      second beta distribution parameter.  q must be gt 0.0.
  51      betai  the incomplete beta function ratio is the probability that a
  52             random variable from a beta distribution having parameters
  53             p and q will be less than or equal to x.
  54   */
  55   double c, finsum, p, ps, q, term, xb, xi, y, dbetai, p1;
  56   static double eps = 0.0, alneps = 0.0, sml = 0.0, alnsml = 0.0;
  57   int i, n, ib;
  58
  59   /* I'm not sure these tolerances are optimal */
  60   if (eps == 0.0) {
  61     eps = macheps();
  62     alneps = log(eps);
  63     sml = macheps();
  64     alnsml = log(sml);
  65   }
  66
  67   y = x;
  68   p = pin;
  69   q = qin;
  70   if (q > p || x >= 0.8)
  71     if (x >= 0.2) {
  72       y = 1.0 - y;
  73       p = qin;
  74       q = pin;
  75     }
  76
  77   if ((p + q) * y / (p + 1.0) < eps) {
  78     dbetai = 0.0;
  79     xb = p * log(Max(y, sml)) - log(p) - logbeta(p, q);
  80     if (xb > alnsml && y != 0.0) dbetai = exp(xb);
  81     if (y != x || p != pin) dbetai = 1.0 - dbetai;
  82   }
  83   else {
  84     /*
  85      *  evaluate the infinite sum first.  term will equal
  86      *  y**p/beta(ps,p) * (1.-ps)-sub-i * y**i / fac(i) .
  87      */
  88     ps = q - floor(q);
  89     if (ps == 0.0) ps = 1.0;
  90     xb = p * log(y) - logbeta(ps, p) - log(p);
  91     dbetai = 0.0;
  92     if (xb >= alnsml) {
  93
  94       dbetai = exp(xb);
  95       term = dbetai * p;
  96       if (ps != 1.0) {
  97         n = Max(alneps / log(y), 4.0);
  98         for (i = 1; i <= n; i++) {
  99           xi = i;
 100           term = term * (xi - ps) * y / xi;
 101           dbetai = dbetai + term / (p + xi);
 102         }
 103       }
 104     }
 105     /*
 106      * now evaluate the finite sum, maybe.
 107      */
 108     if (q > 1.0) {
 109
 110       xb = p * log(y) + q * log(1.0 - y) - logbeta(p,q) - log(q);
 111       ib = Max(xb / alnsml, 0.0);
 112       term = exp(xb - ((float) ib) * alnsml);
 113       c = 1.0 / (1.0 - y);
 114       p1 = q * c / (p + q - 1.0);
 115
 116       finsum = 0.0;
 117       n = q;
 118       if (q == (float) n) n = n - 1;
 119       for (i = 1; i <= n; i++) {
 120         if (p1 <= 1.0 && term / eps <= finsum) break;
 121         xi = i;
 122         term = (q - xi + 1.0) * c * term / (p + q - xi);
 123
 124         if (term > 1.0) ib = ib - 1;
 125         if (term > 1.0) term = term * sml;
 126
 127         if (ib==0) finsum = finsum + term;
 128       }
 129
 130       dbetai = dbetai + finsum;
 131     }
 132     if (y != x || p != pin) dbetai = 1.0 - dbetai;
 133     dbetai = Max(Min(dbetai, 1.0), 0.0);
 134   }
 135   return(dbetai);
 136 }
 137
 138 /*
 139   xinbta.f -- translated by f2c and modified
 140
 141   algorithm as 109 appl. statist. (1977), vol.26, no.1
 142   (replacing algorithm as 64  appl. statist. (1973), vol.22, no.3)
 143
 144   Remark AS R83 has been incorporated in this version.
 145
 146   Computes inverse of the incomplete beta function
 147   ratio for given positive values of the arguments
 148   p and q, alpha between zero and one.
 149   log of complete beta function, beta, is assumed to be known.
 150
 151   Auxiliary function required: betai
 152
 153   SAE below is the most negative decimal exponent which does not
 154   cause an underflow; a value of -308 or thereabouts will often be
 155 */
 156
 157 static double xinbta(p, q, beta, alpha, ifault)
 158      double *p, *q, *beta, *alpha;
 159      int *ifault;
 160 {
 161   /* Initialized data */
 162   static double sae = -30.0; /* this should be sufficient */
 163   static double zero = 0.0;
 164   static double one = 1.0;
 165   static double two = 2.0;
 166   static double three = 3.0;
 167   static double four = 4.0;
 168   static double five = 5.0;
 169   static double six = 6.0;
 170
 171   /* System generated locals */
 172   double ret_val, d_1, d_2;
 173
 174   /* Local variables */
 175   static int indx;
 176   static double prev, a, g, h, r, s, t, w, y, yprev, pp, qq;
 177   static double sq, tx, adj, acu;
 178   static int iex;
 179   static double fpu, xin;
 180
 181   /* Define accuracy and initialise. */
 182   fpu = sae * 10.;
 183   ret_val = *alpha;
 184
 185   /* test for admissibility of parameters */
 186   *ifault = 1;
 187   if (*p <= zero || *q <= zero) return ret_val;
 188   *ifault = 2;
 189   if (*alpha < zero || *alpha > one) return ret_val;
 190   *ifault = 0;
 191   if (*alpha == zero || *alpha == one) return ret_val;
 192
 193   /* change tail if necessary */
 194   if (*alpha <= .5) {
 195     a = *alpha;
 196     pp = *p;
 197     qq = *q;
 198     indx = FALSE;
 199   }
 200   else {
 201     a = one - *alpha;
 202     pp = *q;
 203     qq = *p;
 204     indx = TRUE;
 205   }
 206
 207   /* calculate the initial approximation */
 208   r = sqrt(-log(a * a));
 209   y = r - (r * .27061 + 2.30753) / (one + (r * .04481 + .99229) * r);
 210   if (pp > one && qq > one) {
 211     r = (y * y - three) / six;
 212     s = one / (pp + pp - one);
 213     t = one / (qq + qq - one);
 214     h = two / (s + t);
 215     d_1 = y * sqrt(h + r) / h;
 216     d_2 = (t - s) * (r + five / six - two / (three * h));
 217     w = d_1 - d_2;
 218     ret_val = pp / (pp + qq * exp(w + w));
 219   }
 220   else {
 221     r = qq + qq;
 222     t = one / (qq * 9.);
 223     /* Computing 3rd power */
 224     d_1 = one - t + y * sqrt(t);
 225     t = r * (d_1 * d_1 * d_1);
 226     if (t <= zero) {
 227       ret_val = one - exp((log((one - a) * qq) + *beta) / qq);
 228     }
 229     else {
 230       t = (four * pp + r - two) / t;
 231       if (t <= one) ret_val = exp((log(a * pp) + *beta) / pp);
 232       else ret_val = one - two / (t + one);
 233     }
 234   }
 235
 236
 237   /*
 238     solve for x by a modified newton-raphson method, using the function betai
 239   */
 240   r = one - pp;
 241   t = one - qq;
 242   yprev = zero;
 243   sq = one;
 244   prev = one;
 245   if (ret_val < 1e-4) ret_val = 1e-4;
 246   if (ret_val > .9999) ret_val = .9999;
 247   /* Computing MAX, two 2nd powers */
 248   d_1 = -5.0 / (pp * pp) - 1.0 / (a * a) - 13.0;
 249   iex = (sae > d_1) ? sae : d_1;
 250   acu = pow(10.0, (double) iex);
 251   do {
 252     y = betai(ret_val, pp, qq);
 253     if (*ifault != 0) {
 254       *ifault = 3;
 255       return ret_val;
 256     }
 257     xin = ret_val;
 258     y = (y - a) * exp(*beta + r * log(xin) + t * log(one - xin));
 259     if (y * yprev <= zero) {
 260       prev = (sq > fpu) ? sq : fpu;
 261     }
 262     g = one;
 263     do {
 264       adj = g * y;
 265       sq = adj * adj;
 266       if (sq < prev) {
 267         tx = ret_val - adj;
 268         if (tx >= zero && tx <= one) {
 269           if (prev <= acu || y * y <= acu) {
 270             if (indx) ret_val = one - ret_val;
 271             return ret_val;
 272           }
 273           if (tx != zero && tx != one) break;
 274         }
 275       }
 276       g /= three;
 277     } while (TRUE);
 278     if (tx == ret_val) {
 279       if (indx) ret_val = one - ret_val;
 280       return ret_val;
 281     }
 282     ret_val = tx;
 283     yprev = y;
 284   } while (TRUE);
 285   return ret_val;
 286 }