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