lib/cdists.c

   1 #include "xmath.h"
   2
   3 extern void normbase(double *,double *);
   4 extern void gammabase(double *,double *, double *);
   5 extern void studentbase(double *,double *, double *);
   6 extern void betabase(double *, double *,double *,int *, int *,double *);
   7 extern int max(int, int);
   8 extern int min(int, int);
   9
  10 extern void xlfail(char *);
  11
  12 extern double ppnd(), gamma(), bivnor(), uni(), ppgamma(), ppbeta(),
  13   ppstudent();
  14
  15 /* forward declaration */
  16 extern double tdens();
  17
  18 #ifndef PI
  19 #define PI 3.14159265358979323846
  20 #endif /* PI*/
  21 #define TRUE 1
  22 #define FALSE 0
  23
  24 /*
  25  * Under ULTRIX 3.1 (the cc1.31 compilers in particular) the _G0 math
  26  * library does not really exist! You either need to figure out how
  27  * to get tan() and floor() in at load time for kcl, or use the ones
  28  * here.
  29  */
  30
  31 #ifdef mips
  32 #ifdef UX31
  33 double tan(x)
  34      double x;
  35 {
  36   return(sin(x) / cos(x));
  37 }
  38 double floor(x)
  39      double x;
  40 {
  41   long ix = x;
  42   double dx = ix;
  43   return((dx <= x) ? dx : dx - 1.0);
  44 }
  45 #endif
  46 #endif
  47
  48 static void
  49 checkflag(int flag)
  50 {
  51   /* do nothing for now */
  52 }
  53
  54 static void
  55 checkexp(double a)
  56 {
  57   if (a <= 0.0) {
  58     xlfail("non-positive gamma or beta exponent");
  59   }
  60 }
  61
  62 static void
  63 checkdf(double df)
  64 {
  65   if (df <= 0.0) {
  66     xlfail("non-positive degrees of freedom");
  67   }
  68 }
  69
  70 static void
  71 checkprob(double p, int zerostrict, int onestrict)
  72 {
  73   if (zerostrict) {
  74     if (p <= 0.0) xlfail("non-positive probability argument");
  75   } else {
  76     if (p < 0.0) xlfail("negative probability argument");
  77   }
  78   if (onestrict) {
  79     if (p >= 1.0) xlfail("probability argument not less than one");
  80   } else {
  81     if (p > 1.0) xlfail("probability argument greater than one");
  82   }
  83 }
  84
  85 static void
  86 checkrho(double r)
  87 {
  88   if (r < -1 || r > 1) {
  89     xlfail("correlation out of range");
  90   }
  91 }
  92
  93 static void
  94 checkpoisson(double L)
  95 {
  96   if (L < 0.0) {
  97     xlfail("negative Poisson mean");
  98   }
  99 }
 100
 101 static void
 102 checkbinomial(int n, double p)
 103 {
 104   if (p < 0.0 || p > 1.0) xlfail("binomial p out of range");
 105   if (n < 1) xlfail("non-positive binomial n");
 106 }
 107
 108 /*****************************************************************************/
 109 /*****************************************************************************/
 110 /**                                                                         **/
 111 /**                          Uniform Distribution                           **/
 112 /**                                                                         **/
 113 /*****************************************************************************/
 114 /*****************************************************************************/
 115
 116 /* uniform generator - avoids zero and one */
 117 double unirand()
 118 {
 119   double u;
 120   do {
 121     u = uni();
 122   } while ((u <= 0.0) || (u >= 1.0));
 123   return(u);
 124 }
 125
 126 /*****************************************************************************/
 127 /*****************************************************************************/
 128 /**                                                                         **/
 129 /**                          Normal Distribution                            **/
 130 /**                                                                         **/
 131 /*****************************************************************************/
 132 /*****************************************************************************/
 133
 134 /* standard normal cdf */
 135 double normalcdf(x)
 136      double x;
 137 {
 138   double p;
 139   normbase(&x, &p);
 140   return(p);
 141 }
 142
 143 /* standard normal quantile function */
 144 double normalquant(p)
 145      double p;
 146 {
 147   int flag;
 148   double x;
 149
 150   checkprob(p, TRUE, TRUE);
 151   x = ppnd(p, &flag);
 152   checkflag(flag);
 153   return(x);
 154 }
 155
 156 /* standard normal density */
 157 double normaldens(x)
 158      double x;
 159 {
 160   return(exp(- 0.5 * x * x) / sqrt(2.0 * PI));
 161 }
 162
 163 /* standard normal generator */
 164 double normalrand()
 165 {
 166   double x, y, u, u1, v;
 167   static double c = -1.0;
 168
 169   if (c < 0.0) c = sqrt(2.0 / exp(1.0));
 170
 171   /* ratio of uniforms with linear pretest */
 172   do {
 173     u = unirand();
 174     u1 = unirand();
 175     v = c * (2 * u1 - 1);
 176     x = v / u;
 177     y = x * x / 4.0;
 178   } while(y > (1 - u) && y > - log(u));
 179   return(x);
 180 }
 181
 182 /* bivariate normal cdf */
 183 double bnormcdf(x, y, r)
 184   double x, y, r;
 185 {
 186   checkrho(r);
 187   return(bivnor(-x, -y, r));
 188 }
 189
 190 /*****************************************************************************/
 191 /*****************************************************************************/
 192 /**                                                                         **/
 193 /**                          Cauchy Distribution                            **/
 194 /**                                                                         **/
 195 /*****************************************************************************/
 196 /*****************************************************************************/
 197
 198 double
 199 cauchycdf(double dx)
 200 {
 201   return((atan(dx) + PI / 2) / PI);
 202 }
 203
 204 double
 205 cauchyquant(double p)
 206 {
 207   checkprob(p, TRUE, TRUE);
 208   return(tan(PI * (p - 0.5)));
 209 }
 210
 211 double
 212 cauchydens(double dx)
 213 {
 214   return(tdens(dx, 1.0));
 215 }
 216
 217 /* cauchy generator */
 218 double
 219 cauchyrand()
 220 {
 221   double u1, u2, v1, v2;
 222
 223   /* ratio of uniforms on half disk */
 224   do {
 225     u1 = unirand();
 226     u2 = unirand();
 227     v1 = 2.0 * u1 - 1.0;
 228     v2 = u2;
 229   } while(v1 * v1 + v2 * v2 > 1.0);
 230   return(v1 / v2);
 231 }
 232
 233 /*****************************************************************************/
 234 /**                                                                         **/
 235 /**                           Gamma Distribution                            **/
 236 /**                                                                         **/
 237 /*****************************************************************************/
 238
 239 double
 240 gammacdf(double x, double a)
 241 {
 242   double p;
 243
 244   checkexp(a);
 245   if (x <= 0.0) p = 0.0;
 246   else gammabase(&x, &a, &p);
 247   return(p);
 248 }
 249
 250 double
 251 gammaquant(double p, double a)
 252 {
 253   int flag;
 254   double x;
 255
 256   checkexp(a);
 257   checkprob(p, FALSE, TRUE);
 258   x = ppgamma(p, a, &flag);
 259   checkflag(flag);
 260   return(x);
 261 }
 262
 263 /* gamma density */
 264 double gammadens(double x, double a)
 265 {
 266   double dens;
 267
 268   checkexp(a);
 269   if (x <= 0.0) {
 270     dens = 0.0;
 271   } else {
 272     dens = exp(log(x) * (a - 1) - x - gamma(a));
 273   }
 274   return(dens);
 275 }
 276
 277 /* gamma generator */
 278 double
 279 gammarand(double a)
 280 {
 281   double x, u0, u1, u2, v, w, c, c1, c2, c3, c4, c5;
 282   static double e = -1.0;
 283   int done;
 284
 285   checkexp(a);
 286   if (e < 0.0) e = exp(1.0);
 287
 288   if (a < 1.0) {
 289     /* Ahrens and Dieter algorithm */
 290     done = FALSE;
 291     c = (a + e) / e;
 292     do {
 293       u0 = unirand();
 294       u1 = unirand();
 295       v = c * u0;
 296       if (v <= 1.0) {
 297         x = exp(log(v) / a);
 298         if (u1 <= exp(-x)) done = TRUE;
 299       }
 300       else {
 301         x = -log((c - v) / a);
 302         if (x > 0.0 && u1 < exp((a - 1.0) * log(x))) done = TRUE;
 303       }
 304     } while(! done);
 305   }
 306   else if (a == 1.0) x = -log(unirand());
 307   else {
 308     /* Cheng and Feast algorithm */
 309     c1 = a - 1.0;
 310     c2 = (a - 1.0 / (6.0 * a)) / c1;
 311     c3 = 2.0 / c1;
 312     c4 = 2.0 / (a - 1.0) + 2.0;
 313     c5 = 1.0 / sqrt(a);
 314     do {
 315       do {
 316         u1 = unirand();
 317         u2 = unirand();
 318         if (a > 2.5) u1 = u2 + c5 * (1.0 - 1.86 * u1);
 319       } while (u1 <= 0.0 || u1 >= 1.0);
 320       w = c2 * u2 / u1;
 321     } while ((c3 * u1 + w + 1.0/w) > c4 && (c3 * log(u1) - log(w) + w) > 1.0);
 322     x = c1 * w;
 323   }
 324   return(x);
 325 }
 326
 327 /*****************************************************************************/
 328 /**                                                                         **/
 329 /**                        Chi-Square Distribution                          **/
 330 /**                                                                         **/
 331 /*****************************************************************************/
 332
 333 double chisqcdf(double x, double df)
 334 {
 335   double p, a;
 336
 337   checkdf(df);
 338   a = 0.5 * df;
 339   x = 0.5 * x;
 340   if (x <= 0.0) {
 341     p = 0.0;
 342   } else {
 343     gammabase(&x, &a, &p);
 344   }
 345   return(p);
 346 }
 347
 348 double
 349 chisqquant(double p, double df)
 350 {
 351   double x, a;
 352   int flag;
 353
 354   checkdf(df);
 355   checkprob(p, FALSE, TRUE);
 356   a = 0.5 * df;
 357   x = 2.0 * ppgamma(p, a, &flag);
 358   checkflag(flag);
 359   return(x);
 360 }
 361
 362 double
 363 chisqdens(double dx, double da)
 364 {
 365   checkdf(da);
 366   da = 0.5 * da;
 367   dx = 0.5 * dx;
 368   return(0.5 * gammadens(dx, da));
 369 }
 370
 371 double
 372 chisqrand(double df)
 373 {
 374   checkdf(df);
 375   return(2.0 * gammarand(df / 2.0));
 376 }
 377
 378 /*****************************************************************************/
 379 /**                                                                         **/
 380 /**                           Beta Distribution                             **/
 381 /**                                                                         **/
 382 /*****************************************************************************/
 383
 384 double
 385 betacdf(double x, double a, double b)
 386 {
 387   double p;
 388   int ia, ib;
 389
 390   checkexp(a); checkexp(b);
 391   ia = a; ib = b;
 392   if (x <= 0.0) {
 393     p = 0.0;
 394   } else {
 395     if (x >= 1.0) {
 396       p = 1.0;
 397     } else {
 398       betabase(&x, &a, &b, &ia, &ib, &p);
 399     }
 400   }
 401   return(p);
 402 }
 403
 404 double
 405 betaquant(double p, double a, double b)
 406 {
 407   double x;
 408   int flag;
 409
 410   checkexp(a);
 411   checkexp(b);
 412   checkprob(p, FALSE, FALSE);
 413   x = ppbeta(p, a, b, &flag);
 414   checkflag(flag);
 415   return(x);
 416 }
 417
 418 static double
 419 logbeta(double a, double b)
 420 {
 421   static double da = 0.0, db = 0.0, lbeta = 0.0;
 422
 423   if (a != da || b != db) { /* cache most recent call */
 424     da = a; db = b;
 425     lbeta = gamma(da) + gamma(db) - gamma(da + db);
 426   }
 427   return(lbeta);
 428 }
 429
 430 /* beta density */
 431 double
 432 betadens(double x, double a, double b)
 433 {
 434   double dens;
 435
 436   checkexp(a);
 437   checkexp(b);
 438   if (x <= 0.0 || x >= 1.0) {
 439     dens = 0.0;
 440   } else {
 441     dens = exp(log(x) * (a - 1) + log(1 - x) * (b - 1) - logbeta(a, b));
 442   }
 443   return(dens);
 444 }
 445
 446 /* beta generator */
 447 double
 448 betarand(double a, double b)
 449 {
 450   double x, y;
 451
 452   checkexp(a);
 453   checkexp(b);
 454   x = gammarand(a);
 455   y = gammarand(b);
 456   return(x / (x + y));
 457 }
 458
 459 /*****************************************************************************/
 460 /**                                                                         **/
 461 /**                            t Distribution                               **/
 462 /**                                                                         **/
 463 /*****************************************************************************/
 464
 465 /* t cdf */
 466 double
 467 tcdf(double x, double df)
 468 {
 469   double p;
 470
 471   checkdf(df);
 472   studentbase(&x, &df, &p);
 473   return(p);
 474 }
 475
 476 double
 477 tquant(double p, double df)
 478 {
 479   double x;
 480   int flag;
 481
 482   checkdf(df);
 483   checkprob(p, TRUE, TRUE);
 484   x = ppstudent(p, df, &flag);
 485   checkflag(flag);
 486   return(x);
 487 }
 488
 489 double
 490 tdens(double x, double a)
 491 {
 492   double dens;
 493
 494   checkdf(a);
 495   dens = (1.0 / sqrt(a * PI))
 496        * exp(gamma(0.5 * (a + 1)) - gamma(0.5 * a)
 497              - 0.5 * (a + 1) * log(1.0 + x * x / a));
 498   return(dens);
 499 }
 500
 501 double trand(double df)
 502 {
 503   checkdf(df);
 504   return(normalrand() / sqrt(chisqrand(df) / df));
 505 }
 506
 507 /*****************************************************************************/
 508 /**                                                                         **/
 509 /**                            F Distribution                               **/
 510 /**                                                                         **/
 511 /*****************************************************************************/
 512
 513 double fcdf(double x, double ndf, double ddf)
 514 {
 515   double p, a, b;
 516
 517   checkdf(ndf); checkdf(ddf);
 518   a = 0.5 * ddf;
 519   b = 0.5 * ndf;
 520   if (x <= 0.0) {
 521     p = 0.0;
 522   } else {
 523     x = a / (a + b * x);
 524     p = 1.0 - betacdf(x, a, b);
 525   }
 526   return(p);
 527 }
 528
 529 double fquant(double p, double ndf, double  ddf)
 530 {
 531   double x, a, b;
 532   int flag;
 533
 534   checkdf(ndf); checkdf(ddf);
 535   checkprob(p, FALSE, TRUE);
 536   a = 0.5 * ddf;
 537   b = 0.5 * ndf;
 538   if (p == 0.0) {
 539     x = 0.0;
 540   } else {
 541     p = 1.0 - p;
 542     x = ppbeta(p, a, b, &flag);
 543     checkflag(flag);
 544     x = a * (1.0 / x - 1.0) / b;
 545   }
 546   return(x);
 547 }
 548
 549 double
 550 fdens(double dx, double da, double db)
 551 {
 552   double dens;
 553
 554   checkdf(da);
 555   checkdf(db);
 556   if (dx <= 0.0) {
 557     dens = 0.0;
 558   } else {
 559     dens = exp(0.5 * da * log(da) + 0.5 * db *log(db)
 560                + (0.5 * da - 1.0) * log(dx)
 561                - logbeta(0.5 * da, 0.5 * db)
 562                - 0.5 * (da + db) * log(db + da * dx));
 563   }
 564   return(dens);
 565 }
 566
 567 /* f generator */
 568 double frand(double ndf, double ddf)
 569 {
 570   checkdf(ndf);
 571   checkdf(ddf);
 572   return((ddf * chisqrand(ndf)) / (ndf * chisqrand(ddf)));
 573 }
 574
 575 /*****************************************************************************/
 576 /**                                                                         **/
 577 /**                         Poisson Distribution                            **/
 578 /**                                                                         **/
 579 /*****************************************************************************/
 580
 581 static double
 582 poisson_cdf(int k, double L)
 583 {
 584   double dp, dx;
 585
 586   if (k < 0) {
 587     dp = 0.0;
 588   } else {
 589     if (L == 0.0) {
 590       dp = (k < 0) ? 0.0 : 1.0;
 591     } else {
 592       dx = k + 1.0;
 593       gammabase(&L, &dx, &dp);
 594       dp = 1.0 - dp;
 595     }
 596   }
 597   return(dp);
 598 }
 599
 600 double
 601 poissoncdf(double k, double L)
 602 {
 603   checkpoisson(L);
 604   return(poisson_cdf((int) floor(k), L));
 605 }
 606
 607 int
 608 poissonquant(double x, double L)
 609 {
 610   int k, k1, k2, del, ia;
 611   double m, s, p1, p2, pk;
 612
 613   checkpoisson(L);
 614   checkprob(x, FALSE, TRUE);
 615   m = L;
 616   s = sqrt(L);
 617   del = max(1, (int) (0.2 * s));
 618
 619   if (x == 0.0) {
 620     k = 0.0;
 621   } else {
 622     k = m + s * ppnd(x, &ia);
 623   }
 624   k1 = k;
 625   k2 = k;
 626
 627   do {
 628     k1 = k1 - del; k1 = max(0, k1);
 629     p1 = poisson_cdf(k1, L);
 630   } while (k1 > 0 && p1 > x);
 631
 632   if (k1 == 0 && p1 >= x) {
 633     return(k1);
 634   }
 635
 636   do {
 637     k2 = k2 + del;
 638     p2 = poisson_cdf(k2, L);
 639   } while (p2 < x);
 640
 641   while (k2 - k1 > 1) {
 642     k = (k1 + k2) / 2;
 643     pk = poisson_cdf(k, L);
 644     if (pk < x) {
 645       k1 = k; p1 = pk;
 646     } else {
 647       k2 = k; p2 = pk;
 648     }
 649   }
 650   return(k2);
 651 }
 652
 653 double
 654 poissonpmf(int k, double L)
 655 {
 656   double dx, dp;
 657
 658   checkpoisson(L);
 659   dx = k;
 660   if (L == 0.0) {
 661     dp = (k == 0) ? 1.0 : 0.0;
 662   }  else {
 663     if (dx < 0.0) {
 664       dp = 0.0;
 665     } else {
 666       dp = exp(dx * log(L) - L - gamma(dx + 1.0));
 667     }
 668   }
 669   return(dp);
 670 }
 671
 672 /* poisson random generator from Numerical Recipes */
 673 int poissonrand(double xm)
 674 {
 675   static double sqrt2xm, logxm, expxm, g, oldxm = -1.0;
 676   double t, y;
 677   int k;
 678
 679   checkpoisson(xm);
 680   if (xm < 12.0) {
 681     if (xm != oldxm) { expxm = exp(-xm); oldxm = xm; }
 682     k = -1;
 683     t = 1.0;
 684     do {
 685       k++;
 686       t *= uni();
 687     } while (t > expxm);
 688   } else {
 689     if (xm != oldxm) {
 690       oldxm = xm;
 691       sqrt2xm = sqrt(2.0 * xm);
 692       logxm = log(xm);
 693       g = xm * logxm - gamma(xm + 1.0);
 694     }
 695     do {
 696       do {
 697         y = tan(PI * uni());
 698         k = floor(sqrt2xm * y + xm);
 699       } while (k < 0);
 700       t = 0.9 * (1.0 + y * y) * exp(k * logxm - gamma(k + 1.0) - g);
 701     } while (uni() > t);
 702   }
 703   return (k);
 704 }
 705
 706 /*****************************************************************************/
 707 /**                                                                         **/
 708 /**                        Binomial Distribution                            **/
 709 /**                                                                         **/
 710 /*****************************************************************************/
 711
 712 static double
 713 binomial_cdf(int k, int n, double p)
 714 {
 715   double da, db, dp;
 716   int ia, ib;
 717
 718   if (k < 0) {
 719     dp = 0.0;
 720   } else {
 721     if (k >= n) {
 722       dp = 1.0;
 723     } else {
 724       if (p == 0.0) {
 725         dp = (k < 0) ? 0.0 : 1.0;
 726       } else {
 727         if (p == 1.0) {
 728           dp = (k < n) ? 0.0 : 1.0;
 729         } else {
 730           da = k + 1.0;
 731           db = n - k;
 732           ia = floor(da);
 733           ib = floor(db);
 734           betabase(&p, &da, &db, &ia, &ib, &dp);
 735           dp = 1.0 - dp;
 736         }
 737       }
 738     }
 739   }
 740   return(dp);
 741 }
 742
 743 double
 744 binomialcdf(double k, int n, double p)
 745 {
 746   checkbinomial(n, p);
 747   return(binomial_cdf((int) floor(k), n, p));
 748
 749 }
 750
 751 int
 752 binomialquant(double x, int n, double p)
 753 {
 754   int k, k1, k2, del, ia;
 755   double m, s, p1, p2, pk;
 756
 757   checkbinomial(n, p);
 758   checkprob(x, FALSE, FALSE);
 759
 760   m = n * p;
 761   s = sqrt(n * p * (1 - p));
 762   del = max(1, (int) (0.2 * s));
 763
 764   if (x == 0.0) k = 0.0;
 765   else if (x == 1.0) k = n;
 766   else k = m + s * ppnd(x, &ia);
 767   k1 = k; k2 = k;
 768
 769   do {
 770     k1 = k1 - del; k1 = max(0, k1);
 771     p1 = binomial_cdf(k1, n, p);
 772   } while (k1 > 0 && p1 > p);
 773   if (k1 == 0 && p1 >= x) return(k1);
 774
 775   do {
 776     k2 = k2 + del; k2 = min(n, k2);
 777     p2 = binomial_cdf(k2, n, p);
 778   } while (k2 < n && p2 < x);
 779   if (k2 == n && p2 <= x) return(k2);
 780
 781   while (k2 - k1 > 1) {
 782     k = (k1 + k2) / 2;
 783     pk = binomial_cdf(k, n, p);
 784     if (pk < x) { k1 = k; p1 = pk; }
 785     else { k2 = k; p2 = pk; }
 786   }
 787   return(k2);
 788 }
 789
 790 double binomialpmf(k, n, p)
 791      int k, n;
 792      double p;
 793 {
 794   double dx, dp;
 795
 796   checkbinomial(n, p);
 797   dx = k;
 798   if (p == 0.0) dp = (k == 0) ? 1.0 : 0.0;
 799   else if (p == 1.0) dp = (k == n) ? 1.0 : 0.0;
 800   else if (dx < 0.0 || dx > n) dp = 0.0;
 801   else dp = exp(gamma(n + 1.0) - gamma(dx + 1.0) - gamma(n - dx + 1.0)
 802                 + dx * log(p) + (n - dx) * log(1.0 - p));
 803   return(dp);
 804 }
 805
 806 /* binomial random generator from Numerical Recipes */
 807 int binomialrand(n, pp)
 808         int n;
 809         double pp;
 810 {
 811   int j, k;
 812   static int nold = -1;
 813   double am, em, g, p, sq, t, y;
 814   static double pold = -1.0, pc, plog, pclog, en, oldg;
 815
 816   checkbinomial(n, pp);
 817
 818   p = (pp <= 0.5) ? pp : 1.0 - pp;
 819
 820   am = n * p;
 821   if (p == 0.0) k = 0;
 822   else if (p == 1.0) k = n;
 823   else if (n < 50) {
 824     k = 0;
 825     for (j = 0; j < n; j++) if (uni() < p) k++;
 826   }
 827   else if (am < 1.0) {
 828     g = exp(-am);
 829     t = 1.0;
 830     k = -1;
 831     do {
 832       k++;
 833       t *= uni();
 834     } while (t > g);
 835     if (k > n) k = n;
 836   }
 837   else {
 838     if (n != nold) {
 839       en = n;
 840       oldg = gamma(en + 1.0);
 841       nold = n;
 842     }
 843     if (p != pold) {
 844       pc = 1.0 - p;
 845       plog = log(p);
 846       pclog = log(pc);
 847       pold = p;
 848     }
 849     sq = sqrt(2.0 * am * pc);
 850     do {
 851       do {
 852         y = tan(PI * uni());
 853         em = sq * y + am;
 854       } while (em < 0.0 || em >= en + 1.0);
 855       em = floor(em);
 856       t = 1.2 * sq * (1.0 + y * y)
 857         * exp(oldg - gamma(em + 1.0) - gamma(en - em + 1.0)
 858               + em * plog + (en - em) * pclog);
 859     } while (uni() > t);
 860     k = em;
 861   }
 862   if (p != pp) k = n - k;
 863   return(k);
 864 }