lsbasics has some numerical processing, essential structure for XLS.
[CommonLispStat.git] / lib / cdists.c
blob7be43792838736766008f48532d2c43ebcb8d5d5
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 *,int *, int *,double *);
7 extern int max(int, int);
8 extern int min(int, int);
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 #define TRUE 1
22 #define FALSE 0
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.
31 #ifdef mips
32 #ifdef UX31
33 double tan(x)
34 double x;
36 return(sin(x) / cos(x));
38 double floor(x)
39 double x;
41 long ix = x;
42 double dx = ix;
43 return((dx <= x) ? dx : dx - 1.0);
45 #endif
46 #endif
48 static void
49 checkflag(int flag)
51 /* do nothing for now */
54 static void
55 checkexp(double a)
57 if (a <= 0.0) {
58 xlfail("non-positive gamma or beta exponent");
62 static void
63 checkdf(double df)
65 if (df <= 0.0) {
66 xlfail("non-positive degrees of freedom");
70 static void
71 checkprob(double p, int zerostrict, int onestrict)
73 if (zerostrict) {
74 if (p <= 0.0) xlfail("non-positive probability argument");
75 } else {
76 if (p < 0.0) xlfail("negative probability argument");
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");
85 static void
86 checkrho(double r)
88 if (r < -1 || r > 1) {
89 xlfail("correlation out of range");
93 static void
94 checkpoisson(double L)
96 if (L < 0.0) {
97 xlfail("negative Poisson mean");
101 static void
102 checkbinomial(int n, double p)
104 if (p < 0.0 || p > 1.0) xlfail("binomial p out of range");
105 if (n < 1) xlfail("non-positive binomial n");
108 /*****************************************************************************/
109 /*****************************************************************************/
110 /** **/
111 /** Uniform Distribution **/
112 /** **/
113 /*****************************************************************************/
114 /*****************************************************************************/
116 /* uniform generator - avoids zero and one */
117 double unirand()
119 double u;
120 do {
121 u = uni();
122 } while ((u <= 0.0) || (u >= 1.0));
123 return(u);
126 /*****************************************************************************/
127 /*****************************************************************************/
128 /** **/
129 /** Normal Distribution **/
130 /** **/
131 /*****************************************************************************/
132 /*****************************************************************************/
134 /* standard normal cdf */
135 double normalcdf(x)
136 double x;
138 double p;
139 normbase(&x, &p);
140 return(p);
143 /* standard normal quantile function */
144 double normalquant(p)
145 double p;
147 int flag;
148 double x;
150 checkprob(p, TRUE, TRUE);
151 x = ppnd(p, &flag);
152 checkflag(flag);
153 return(x);
156 /* standard normal density */
157 double normaldens(x)
158 double x;
160 return(exp(- 0.5 * x * x) / sqrt(2.0 * PI));
163 /* standard normal generator */
164 double normalrand()
166 double x, y, u, u1, v;
167 static double c = -1.0;
169 if (c < 0.0) c = sqrt(2.0 / exp(1.0));
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);
182 /* bivariate normal cdf */
183 double bnormcdf(x, y, r)
184 double x, y, 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 int 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 int 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(int 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((int) floor(k), L));
608 poissonquant(double x, double L)
610 int 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, (int) (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(int 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 int poissonrand(double xm)
675 static double sqrt2xm, logxm, expxm, g, oldxm = -1.0;
676 double t, y;
677 int 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(int k, int n, double p)
715 double da, db, dp;
716 int 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, int n, double p)
746 checkbinomial(n, p);
747 return(binomial_cdf((int) floor(k), n, p));
752 binomialquant(double x, int n, double p)
754 int k, k1, k2, del, ia;
755 double m, s, p1, p2, pk;
757 checkbinomial(n, p);
758 checkprob(x, FALSE, FALSE);
760 m = n * p;
761 s = sqrt(n * p * (1 - p));
762 del = max(1, (int) (0.2 * s));
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;
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);
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);
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; }
787 return(k2);
790 double binomialpmf(k, n, p)
791 int k, n;
792 double p;
794 double dx, dp;
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);
806 /* binomial random generator from Numerical Recipes */
807 int binomialrand(n, pp)
808 int n;
809 double pp;
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;
816 checkbinomial(n, pp);
818 p = (pp <= 0.5) ? pp : 1.0 - pp;
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++;
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;
837 else {
838 if (n != nold) {
839 en = n;
840 oldg = gamma(en + 1.0);
841 nold = n;
843 if (p != pold) {
844 pc = 1.0 - p;
845 plog = log(p);
846 pclog = log(pc);
847 pold = p;
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;
862 if (p != pp) k = n - k;
863 return(k);