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 *,long *, long *,double *);
   7 extern long max(long, long);
   8 extern long min(long, long);
   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 #ifdef __LP64__
  22 #define TRUE 1L
  23 #define FALSE 0L
  24 #else
  25 #define TRUE 1
  26 #define FALSE 0
  27 #endif
  28
  29 /*
  30  * Under ULTRIX 3.1 (the cc1.31 compilers in particular) the _G0 math
  31  * library does not really exist! You either need to figure out how
  32  * to get tan() and floor() in at load time for kcl, or use the ones
  33  * here.
  34  */
  35
  36 #ifdef mips
  37 #ifdef UX31
  38 double tan(double x)
  39 {
  40   return(sin(x) / cos(x));
  41 }
  42
  43 double floor(double x)
  44 {
  45   long ix = x;
  46   double dx = ix;
  47   return((dx <= x) ? dx : dx - 1.0);
  48 }
  49 #endif
  50 #endif
  51
  52 static void
  53 checkflag(int flag)
  54 {
  55   /* do nothing for now */
  56 }
  57
  58 static void
  59 checkexp(double a)
  60 {
  61   if (a <= 0.0) {
  62     xlfail("non-positive gamma or beta exponent");
  63   }
  64 }
  65
  66 static void
  67 checkdf(double df)
  68 {
  69   if (df <= 0.0) {
  70     xlfail("non-positive degrees of freedom");
  71   }
  72 }
  73
  74 static void
  75 checkprob(double p, long zerostrict, long onestrict)
  76 {
  77   if (zerostrict) {
  78     if (p <= 0.0) xlfail("non-positive probability argument");
  79   } else {
  80     if (p < 0.0) xlfail("negative probability argument");
  81   }
  82   if (onestrict) {
  83     if (p >= 1.0) xlfail("probability argument not less than one");
  84   } else {
  85     if (p > 1.0) xlfail("probability argument greater than one");
  86   }
  87 }
  88
  89 static void
  90 checkrho(double r)
  91 {
  92   if (r < -1 || r > 1) {
  93     xlfail("correlation out of range");
  94   }
  95 }
  96
  97 static void
  98 checkpoisson(double L)
  99 {
 100   if (L < 0.0) {
 101     xlfail("negative Poisson mean");
 102   }
 103 }
 104
 105 static void
 106 checkbinomial(long n, double p)
 107 {
 108   if (p < 0.0 || p > 1.0) xlfail("binomial p out of range");
 109   if (n < 1) xlfail("non-positive binomial n");
 110 }
 111
 112 /*****************************************************************************/
 113 /*****************************************************************************/
 114 /**                                                                         **/
 115 /**                          Uniform Distribution                           **/
 116 /**                                                                         **/
 117 /*****************************************************************************/
 118 /*****************************************************************************/
 119
 120 /* uniform generator - avoids zero and one */
 121 double unirand()
 122 {
 123   double u;
 124   do {
 125     u = uni();
 126   } while ((u <= 0.0) || (u >= 1.0));
 127   return(u);
 128 }
 129
 130 /*****************************************************************************/
 131 /*****************************************************************************/
 132 /**                                                                         **/
 133 /**                          Normal Distribution                            **/
 134 /**                                                                         **/
 135 /*****************************************************************************/
 136 /*****************************************************************************/
 137
 138 /* standard normal cdf */
 139 double normalcdf(double x)
 140 {
 141   double p;
 142   normbase(&x, &p);
 143   return(p);
 144 }
 145
 146 /* standard normal quantile function */
 147 double normalquant(double p)
 148 {
 149   int flag;
 150   double x;
 151
 152   checkprob(p, TRUE, TRUE);
 153   x = ppnd(p, &flag);
 154   checkflag(flag);
 155   return(x);
 156 }
 157
 158 /* standard normal density */
 159 double normaldens(double x)
 160 {
 161   return(exp(- 0.5 * x * x) / sqrt(2.0 * PI));
 162 }
 163
 164 /* standard normal generator */
 165 double normalrand()
 166 {
 167   double x, y, u, u1, v;
 168   static double c = -1.0;
 169
 170   if (c < 0.0) c = sqrt(2.0 / exp(1.0));
 171
 172   /* ratio of uniforms with linear pretest */
 173   do {
 174     u = unirand();
 175     u1 = unirand();
 176     v = c * (2 * u1 - 1);
 177     x = v / u;
 178     y = x * x / 4.0;
 179   } while(y > (1 - u) && y > - log(u));
 180   return(x);
 181 }
 182
 183 /* bivariate normal cdf */
 184 double bnormcdf(double x, double y, double 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   long 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   long 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(long 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((long) floor(k), L));
 605 }
 606
 607 long
 608 poissonquant(double x, double L)
 609 {
 610   long 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, (long) (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(long 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 long poissonrand(double xm)
 674 {
 675   static double sqrt2xm, logxm, expxm, g, oldxm = -1.0;
 676   double t, y;
 677   long 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(long k, long n, double p)
 714 {
 715   double da, db, dp;
 716   long 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, long n, double p)
 745 {
 746   checkbinomial(n, p);
 747   return(binomial_cdf((long) floor(k), n, p));
 748 }
 749
 750 long
 751 binomialquant(double x, long n, double p)
 752 {
 753   long k, k1, k2, del, ia;
 754   double m, s, p1, p2, pk;
 755
 756   checkbinomial(n, p);
 757   checkprob(x, FALSE, FALSE);
 758
 759   m = n * p;
 760   s = sqrt(n * p * (1 - p));
 761   del = max(1, (long) (0.2 * s));
 762
 763   if (x == 0.0) k = 0.0;
 764   else if (x == 1.0) k = n;
 765   else k = m + s * ppnd(x, &ia);
 766   k1 = k; k2 = k;
 767
 768   do {
 769     k1 = k1 - del; k1 = max(0, k1);
 770     p1 = binomial_cdf(k1, n, p);
 771   } while (k1 > 0 && p1 > p);
 772   if (k1 == 0 && p1 >= x) return(k1);
 773
 774   do {
 775     k2 = k2 + del; k2 = min(n, k2);
 776     p2 = binomial_cdf(k2, n, p);
 777   } while (k2 < n && p2 < x);
 778   if (k2 == n && p2 <= x) return(k2);
 779
 780   while (k2 - k1 > 1) {
 781     k = (k1 + k2) / 2;
 782     pk = binomial_cdf(k, n, p);
 783     if (pk < x) { k1 = k; p1 = pk; }
 784     else { k2 = k; p2 = pk; }
 785   }
 786   return(k2);
 787 }
 788
 789 double binomialpmf(long k, long n, double p)
 790 {
 791   double dx, dp;
 792
 793   checkbinomial(n, p);
 794   dx = k;
 795   if (p == 0.0) dp = (k == 0) ? 1.0 : 0.0;
 796   else if (p == 1.0) dp = (k == n) ? 1.0 : 0.0;
 797   else if (dx < 0.0 || dx > n) dp = 0.0;
 798   else dp = exp(gamma(n + 1.0) - gamma(dx + 1.0) - gamma(n - dx + 1.0)
 799                 + dx * log(p) + (n - dx) * log(1.0 - p));
 800   return(dp);
 801 }
 802
 803 /* binomial random generator from Numerical Recipes */
 804 long binomialrand(long n, double pp)
 805 {
 806   long j, k;
 807   static long nold = -1;
 808   double am, em, g, p, sq, t, y;
 809   static double pold = -1.0, pc, plog, pclog, en, oldg;
 810
 811   checkbinomial(n, pp);
 812
 813   p = (pp <= 0.5) ? pp : 1.0 - pp;
 814
 815   am = n * p;
 816   if (p == 0.0) k = 0;
 817   else if (p == 1.0) k = n;
 818   else if (n < 50) {
 819     k = 0;
 820     for (j = 0; j < n; j++) if (uni() < p) k++;
 821   }
 822   else if (am < 1.0) {
 823     g = exp(-am);
 824     t = 1.0;
 825     k = -1;
 826     do {
 827       k++;
 828       t *= uni();
 829     } while (t > g);
 830     if (k > n) k = n;
 831   }
 832   else {
 833     if (n != nold) {
 834       en = n;
 835       oldg = gamma(en + 1.0);
 836       nold = n;
 837     }
 838     if (p != pold) {
 839       pc = 1.0 - p;
 840       plog = log(p);
 841       pclog = log(pc);
 842       pold = p;
 843     }
 844     sq = sqrt(2.0 * am * pc);
 845     do {
 846       do {
 847         y = tan(PI * uni());
 848         em = sq * y + am;
 849       } while (em < 0.0 || em >= en + 1.0);
 850       em = floor(em);
 851       t = 1.2 * sq * (1.0 + y * y)
 852         * exp(oldg - gamma(em + 1.0) - gamma(en - em + 1.0)
 853               + em * plog + (en - em) * pclog);
 854     } while (uni() > t);
 855     k = em;
 856   }
 857   if (p != pp) k = n - k;
 858   return(k);
 859 }