initial work towards getting xarray support into CLS.
[CommonLispStat.git] / lib / cdists.c
blobe8cad371dc1a479d89b2486d70080595c140dc16
1 #include "xmath.h"
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);
10 extern void xlfail(char *);
12 extern double ppnd(), gamma(), bivnor(), uni(), ppgamma(), ppbeta(),
13 ppstudent();
15 /* forward declaration */
16 extern double tdens();
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
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.
36 #ifdef mips
37 #ifdef UX31
38 double tan(double x)
40 return(sin(x) / cos(x));
43 double floor(double x)
45 long ix = x;
46 double dx = ix;
47 return((dx <= x) ? dx : dx - 1.0);
49 #endif
50 #endif
52 static void
53 checkflag(int flag)
55 /* do nothing for now */
58 static void
59 checkexp(double a)
61 if (a <= 0.0) {
62 xlfail("non-positive gamma or beta exponent");
66 static void
67 checkdf(double df)
69 if (df <= 0.0) {
70 xlfail("non-positive degrees of freedom");
74 static void
75 checkprob(double p, long zerostrict, long onestrict)
77 if (zerostrict) {
78 if (p <= 0.0) xlfail("non-positive probability argument");
79 } else {
80 if (p < 0.0) xlfail("negative probability argument");
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");
89 static void
90 checkrho(double r)
92 if (r < -1 || r > 1) {
93 xlfail("correlation out of range");
97 static void
98 checkpoisson(double L)
100 if (L < 0.0) {
101 xlfail("negative Poisson mean");
105 static void
106 checkbinomial(long n, double p)
108 if (p < 0.0 || p > 1.0) xlfail("binomial p out of range");
109 if (n < 1) xlfail("non-positive binomial n");
112 /*****************************************************************************/
113 /*****************************************************************************/
114 /** **/
115 /** Uniform Distribution **/
116 /** **/
117 /*****************************************************************************/
118 /*****************************************************************************/
120 /* uniform generator - avoids zero and one */
121 double unirand()
123 double u;
124 do {
125 u = uni();
126 } while ((u <= 0.0) || (u >= 1.0));
127 return(u);
130 /*****************************************************************************/
131 /*****************************************************************************/
132 /** **/
133 /** Normal Distribution **/
134 /** **/
135 /*****************************************************************************/
136 /*****************************************************************************/
138 /* standard normal cdf */
139 double normalcdf(double x)
141 double p;
142 normbase(&x, &p);
143 return(p);
146 /* standard normal quantile function */
147 double normalquant(double p)
149 int flag;
150 double x;
152 checkprob(p, TRUE, TRUE);
153 x = ppnd(p, &flag);
154 checkflag(flag);
155 return(x);
158 /* standard normal density */
159 double normaldens(double x)
161 return(exp(- 0.5 * x * x) / sqrt(2.0 * PI));
164 /* standard normal generator */
165 double normalrand()
167 double x, y, u, u1, v;
168 static double c = -1.0;
170 if (c < 0.0) c = sqrt(2.0 / exp(1.0));
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);
183 /* bivariate normal cdf */
184 double bnormcdf(double x, double y, double r)
186 checkrho(r);
187 return(bivnor(-x, -y, r));
190 /*****************************************************************************/
191 /*****************************************************************************/
192 /** **/
193 /** Cauchy Distribution **/
194 /** **/
195 /*****************************************************************************/
196 /*****************************************************************************/
198 double
199 cauchycdf(double dx)
201 return((atan(dx) + PI / 2) / PI);
204 double
205 cauchyquant(double p)
207 checkprob(p, TRUE, TRUE);
208 return(tan(PI * (p - 0.5)));
211 double
212 cauchydens(double dx)
214 return(tdens(dx, 1.0));
217 /* cauchy generator */
218 double
219 cauchyrand()
221 double u1, u2, v1, v2;
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);
233 /*****************************************************************************/
234 /** **/
235 /** Gamma Distribution **/
236 /** **/
237 /*****************************************************************************/
239 double
240 gammacdf(double x, double a)
242 double p;
244 checkexp(a);
245 if (x <= 0.0) p = 0.0;
246 else gammabase(&x, &a, &p);
247 return(p);
250 double
251 gammaquant(double p, double a)
253 int flag;
254 double x;
256 checkexp(a);
257 checkprob(p, FALSE, TRUE);
258 x = ppgamma(p, a, &flag);
259 checkflag(flag);
260 return(x);
263 /* gamma density */
264 double gammadens(double x, double a)
266 double dens;
268 checkexp(a);
269 if (x <= 0.0) {
270 dens = 0.0;
271 } else {
272 dens = exp(log(x) * (a - 1) - x - gamma(a));
274 return(dens);
277 /* gamma generator */
278 double
279 gammarand(double a)
281 double x, u0, u1, u2, v, w, c, c1, c2, c3, c4, c5;
282 static double e = -1.0;
283 long done;
285 checkexp(a);
286 if (e < 0.0) e = exp(1.0);
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;
300 else {
301 x = -log((c - v) / a);
302 if (x > 0.0 && u1 < exp((a - 1.0) * log(x))) done = TRUE;
304 } while(! done);
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;
324 return(x);
327 /*****************************************************************************/
328 /** **/
329 /** Chi-Square Distribution **/
330 /** **/
331 /*****************************************************************************/
333 double chisqcdf(double x, double df)
335 double p, a;
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);
345 return(p);
348 double
349 chisqquant(double p, double df)
351 double x, a;
352 int flag;
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);
362 double
363 chisqdens(double dx, double da)
365 checkdf(da);
366 da = 0.5 * da;
367 dx = 0.5 * dx;
368 return(0.5 * gammadens(dx, da));
371 double
372 chisqrand(double df)
374 checkdf(df);
375 return(2.0 * gammarand(df / 2.0));
378 /*****************************************************************************/
379 /** **/
380 /** Beta Distribution **/
381 /** **/
382 /*****************************************************************************/
384 double
385 betacdf(double x, double a, double b)
387 double p;
388 long ia, ib;
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);
401 return(p);
404 double
405 betaquant(double p, double a, double b)
407 double x;
408 int flag;
410 checkexp(a);
411 checkexp(b);
412 checkprob(p, FALSE, FALSE);
413 x = ppbeta(p, a, b, &flag);
414 checkflag(flag);
415 return(x);
418 static double
419 logbeta(double a, double b)
421 static double da = 0.0, db = 0.0, lbeta = 0.0;
423 if (a != da || b != db) { /* cache most recent call */
424 da = a; db = b;
425 lbeta = gamma(da) + gamma(db) - gamma(da + db);
427 return(lbeta);
430 /* beta density */
431 double
432 betadens(double x, double a, double b)
434 double dens;
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));
443 return(dens);
446 /* beta generator */
447 double
448 betarand(double a, double b)
450 double x, y;
452 checkexp(a);
453 checkexp(b);
454 x = gammarand(a);
455 y = gammarand(b);
456 return(x / (x + y));
459 /*****************************************************************************/
460 /** **/
461 /** t Distribution **/
462 /** **/
463 /*****************************************************************************/
465 /* t cdf */
466 double
467 tcdf(double x, double df)
469 double p;
471 checkdf(df);
472 studentbase(&x, &df, &p);
473 return(p);
476 double
477 tquant(double p, double df)
479 double x;
480 int flag;
482 checkdf(df);
483 checkprob(p, TRUE, TRUE);
484 x = ppstudent(p, df, &flag);
485 checkflag(flag);
486 return(x);
489 double
490 tdens(double x, double a)
492 double dens;
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);
501 double trand(double df)
503 checkdf(df);
504 return(normalrand() / sqrt(chisqrand(df) / df));
507 /*****************************************************************************/
508 /** **/
509 /** F Distribution **/
510 /** **/
511 /*****************************************************************************/
513 double fcdf(double x, double ndf, double ddf)
515 double p, a, b;
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);
526 return(p);
529 double fquant(double p, double ndf, double ddf)
531 double x, a, b;
532 int flag;
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;
546 return(x);
549 double
550 fdens(double dx, double da, double db)
552 double dens;
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));
564 return(dens);
567 /* f generator */
568 double frand(double ndf, double ddf)
570 checkdf(ndf);
571 checkdf(ddf);
572 return((ddf * chisqrand(ndf)) / (ndf * chisqrand(ddf)));
575 /*****************************************************************************/
576 /** **/
577 /** Poisson Distribution **/
578 /** **/
579 /*****************************************************************************/
581 static double
582 poisson_cdf(long k, double L)
584 double dp, dx;
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;
597 return(dp);
600 double
601 poissoncdf(double k, double L)
603 checkpoisson(L);
604 return(poisson_cdf((long) floor(k), L));
607 long
608 poissonquant(double x, double L)
610 long k, k1, k2, del, ia;
611 double m, s, p1, p2, pk;
613 checkpoisson(L);
614 checkprob(x, FALSE, TRUE);
615 m = L;
616 s = sqrt(L);
617 del = max(1, (long) (0.2 * s));
619 if (x == 0.0) {
620 k = 0.0;
621 } else {
622 k = m + s * ppnd(x, &ia);
624 k1 = k;
625 k2 = k;
627 do {
628 k1 = k1 - del; k1 = max(0, k1);
629 p1 = poisson_cdf(k1, L);
630 } while (k1 > 0 && p1 > x);
632 if (k1 == 0 && p1 >= x) {
633 return(k1);
636 do {
637 k2 = k2 + del;
638 p2 = poisson_cdf(k2, L);
639 } while (p2 < x);
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;
650 return(k2);
653 double
654 poissonpmf(long k, double L)
656 double dx, dp;
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));
669 return(dp);
672 /* poisson random generator from Numerical Recipes */
673 long poissonrand(double xm)
675 static double sqrt2xm, logxm, expxm, g, oldxm = -1.0;
676 double t, y;
677 long k;
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);
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);
703 return (k);
706 /*****************************************************************************/
707 /** **/
708 /** Binomial Distribution **/
709 /** **/
710 /*****************************************************************************/
712 static double
713 binomial_cdf(long k, long n, double p)
715 double da, db, dp;
716 long ia, ib;
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;
740 return(dp);
743 double
744 binomialcdf(double k, long n, double p)
746 checkbinomial(n, p);
747 return(binomial_cdf((long) floor(k), n, p));
750 long
751 binomialquant(double x, long n, double p)
753 long k, k1, k2, del, ia;
754 double m, s, p1, p2, pk;
756 checkbinomial(n, p);
757 checkprob(x, FALSE, FALSE);
759 m = n * p;
760 s = sqrt(n * p * (1 - p));
761 del = max(1, (long) (0.2 * s));
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;
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);
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);
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; }
786 return(k2);
789 double binomialpmf(long k, long n, double p)
791 double dx, dp;
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);
803 /* binomial random generator from Numerical Recipes */
804 long binomialrand(long n, double pp)
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;
811 checkbinomial(n, pp);
813 p = (pp <= 0.5) ? pp : 1.0 - pp;
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++;
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;
832 else {
833 if (n != nold) {
834 en = n;
835 oldg = gamma(en + 1.0);
836 nold = n;
838 if (p != pold) {
839 pc = 1.0 - p;
840 plog = log(p);
841 pclog = log(pc);
842 pold = p;
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;
857 if (p != pp) k = n - k;
858 return(k);