lib/cdists.c

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