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