3 extern double ppnd(), gamma(), bivnor(), uni(), ppgamma(), ppbeta(),
6 /* forward declaration */
10 #define PI 3.14159265358979323846
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
27 return(sin(x
) / cos(x
));
34 return((dx
<= x
) ? dx
: dx
- 1.0);
39 static checkflag(flag
)
42 /* do nothing for now */
48 if (a
<= 0.0) xlfail("non-positive gamma or beta exponent");
54 if (df
<= 0.0) xlfail("non-positive degrees of freedom");
57 static checkprob(p
, zerostrict
, onestrict
)
59 int zerostrict
, onestrict
;
62 if (p
<= 0.0) xlfail("non-positive probability argument");
65 if (p
< 0.0) xlfail("negative probability argument");
68 if (p
>= 1.0) xlfail("probability argument not less than one");
71 if (p
> 1.0) xlfail("probability argument greater than one");
78 if (r
< -1 || r
> 1) xlfail("correlation out of range");
81 static checkpoisson(L
)
84 if (L
< 0.0) xlfail("negative Poisson mean");
87 static checkbinomial(n
, p
)
91 if (p
< 0.0 || p
> 1.0) xlfail("binomial p out of range");
92 if (n
< 1) xlfail("non-positive binomial n");
95 /*****************************************************************************/
96 /*****************************************************************************/
98 /** Uniform Distribution **/
100 /*****************************************************************************/
101 /*****************************************************************************/
103 /* uniform generator - avoids zero and one */
109 } while ((u
<= 0.0) || (u
>= 1.0));
113 /*****************************************************************************/
114 /*****************************************************************************/
116 /** Normal Distribution **/
118 /*****************************************************************************/
119 /*****************************************************************************/
121 /* standard normal cdf */
130 /* standard normal quantile function */
131 double normalquant(p
)
137 checkprob(p
, TRUE
, TRUE
);
143 /* standard normal density */
147 return(exp(- 0.5 * x
* x
) / sqrt(2.0 * PI
));
150 /* standard normal generator */
153 double x
, y
, u
, u1
, v
;
154 static double c
= -1.0;
156 if (c
< 0.0) c
= sqrt(2.0 / exp(1.0));
158 /* ratio of uniforms with linear pretest */
162 v
= c
* (2 * u1
- 1);
165 } while(y
> (1 - u
) && y
> - log(u
));
169 /* bivariate normal cdf */
170 double bnormcdf(x
, y
, r
)
174 return(bivnor(-x
, -y
, r
));
177 /*****************************************************************************/
178 /*****************************************************************************/
180 /** Cauchy Distribution **/
182 /*****************************************************************************/
183 /*****************************************************************************/
189 return((atan(dx
) + PI
/ 2) / PI
);
192 /* Cauchy quantile function */
193 double cauchyquant(p
)
196 checkprob(p
, TRUE
, TRUE
);
197 return(tan(PI
* (p
- 0.5)));
201 double cauchydens(dx
)
204 return(tdens(dx
, 1.0));
207 /* cauchy generator */
210 double u1
, u2
, v1
, v2
;
212 /* ratio of uniforms on half disk */
218 } while(v1
* v1
+ v2
* v2
> 1.0);
222 /*****************************************************************************/
223 /*****************************************************************************/
225 /** Gamma Distribution **/
227 /*****************************************************************************/
228 /*****************************************************************************/
231 double gammacdf(x
, a
)
237 if (x
<= 0.0) p
= 0.0;
238 else gammabase(&x
, &a
, &p
);
242 double gammaquant(p
, a
)
249 checkprob(p
, FALSE
, TRUE
);
250 x
= ppgamma(p
, a
, &flag
);
256 double gammadens(x
, a
)
262 if (x
<= 0.0) dens
= 0.0;
263 else dens
= exp(log(x
) * (a
- 1) - x
- gamma(a
));
267 /* gamma generator */
271 double x
, u0
, u1
, u2
, v
, w
, c
, c1
, c2
, c3
, c4
, c5
;
272 static double e
= -1.0;
276 if (e
< 0.0) e
= exp(1.0);
279 /* Ahrens and Dieter algorithm */
288 if (u1
<= exp(-x
)) done
= TRUE
;
291 x
= -log((c
- v
) / a
);
292 if (x
> 0.0 && u1
< exp((a
- 1.0) * log(x
))) done
= TRUE
;
296 else if (a
== 1.0) x
= -log(unirand());
298 /* Cheng and Feast algorithm */
300 c2
= (a
- 1.0 / (6.0 * a
)) / c1
;
302 c4
= 2.0 / (a
- 1.0) + 2.0;
308 if (a
> 2.5) u1
= u2
+ c5
* (1.0 - 1.86 * u1
);
309 } while (u1
<= 0.0 || u1
>= 1.0);
311 } while ((c3
* u1
+ w
+ 1.0/w
) > c4
&& (c3
* log(u1
) - log(w
) + w
) > 1.0);
317 /*****************************************************************************/
318 /*****************************************************************************/
320 /** Chi-Square Distribution **/
322 /*****************************************************************************/
323 /*****************************************************************************/
325 double chisqcdf(x
, df
)
331 a
= 0.5 * df
; x
= 0.5 * x
;
332 if (x
<= 0.0) p
= 0.0;
333 else gammabase(&x
, &a
, &p
);
337 double chisqquant(p
, df
)
344 checkprob(p
, FALSE
, TRUE
);
346 x
= 2.0 * ppgamma(p
, a
, &flag
);
351 /* chi-square density */
352 double chisqdens(dx
, da
)
358 return(0.5 * gammadens(dx
, da
));
361 /* chi-square generator */
366 return(2.0 * gammarand(df
/ 2.0));
369 /*****************************************************************************/
370 /*****************************************************************************/
372 /** Beta Distribution **/
374 /*****************************************************************************/
375 /*****************************************************************************/
377 double betacdf(x
, a
, b
)
383 checkexp(a
); checkexp(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
);
391 double betaquant(p
, a
, b
)
397 checkexp(a
); checkexp(b
);
398 checkprob(p
, FALSE
, FALSE
);
399 x
= ppbeta(p
, a
, b
, &flag
);
404 static double logbeta(a
, b
)
407 static double da
= 0.0, db
= 0.0, lbeta
= 0.0;
409 if (a
!= da
|| b
!= db
) { /* cache most recent call */
411 lbeta
= gamma(da
) + gamma(db
) - gamma(da
+ db
);
417 double betadens(x
, a
, 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
));
430 double betarand(a
, b
)
442 /*****************************************************************************/
443 /*****************************************************************************/
445 /** t Distribution **/
447 /*****************************************************************************/
448 /*****************************************************************************/
457 studentbase(&x
, &df
, &p
);
461 /* t quantile function */
469 checkprob(p
, TRUE
, TRUE
);
470 x
= ppstudent(p
, df
, &flag
);
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
));
493 return(normalrand() / sqrt(chisqrand(df
) / df
));
496 /*****************************************************************************/
497 /*****************************************************************************/
499 /** F Distribution **/
501 /*****************************************************************************/
502 /*****************************************************************************/
505 double fcdf(x
, ndf
, ddf
)
510 checkdf(ndf
); checkdf(ddf
);
513 if (x
<= 0.0) p
= 0.0;
516 p
= 1.0 - betacdf(x
, a
, b
);
521 /* f quantile function */
522 double fquant(p
, ndf
, ddf
)
528 checkdf(ndf
); checkdf(ddf
);
529 checkprob(p
, FALSE
, TRUE
);
532 if (p
== 0.0) x
= 0.0;
535 x
= ppbeta(p
, a
, b
, &flag
);
537 x
= a
* (1.0 / x
- 1.0) / b
;
543 double fdens(dx
, da
, 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
));
559 double frand(ndf
, ddf
)
564 return((ddf
* chisqrand(ndf
)) / (ndf
* chisqrand(ddf
)));
567 /*****************************************************************************/
568 /*****************************************************************************/
570 /** Poisson Distribution **/
572 /*****************************************************************************/
573 /*****************************************************************************/
575 static double poisson_cdf(k
, L
)
582 else if (L
== 0.0) dp
= (k
< 0) ? 0.0 : 1.0;
585 gammabase(&L
, &dx
, &dp
);
591 double poissoncdf(k
, L
)
596 return(poisson_cdf((int) floor(k
), L
));
599 int poissonquant(x
, L
)
602 int k
, k1
, k2
, del
, ia
;
603 double m
, s
, p1
, p2
, pk
;
606 checkprob(x
, FALSE
, TRUE
);
609 del
= max(1, (int) (0.2 * s
));
611 if (x
== 0.0) k
= 0.0;
612 else k
= m
+ s
* ppnd(x
, &ia
);
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
);
623 p2
= poisson_cdf(k2
, L
);
626 while (k2
- k1
> 1) {
628 pk
= poisson_cdf(k
, L
);
629 if (pk
< x
) { k1
= k
; p1
= pk
; }
630 else { k2
= k
; p2
= pk
; }
635 double poissonpmf(k
, L
)
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));
649 /* poisson random generator from Numerical Recipes */
653 static double sqrt2xm
, logxm
, expxm
, g
, oldxm
= -1.0;
659 if (xm
!= oldxm
) { expxm
= exp(-xm
); oldxm
= xm
; }
670 sqrt2xm
= sqrt(2.0 * xm
);
672 g
= xm
* logxm
- gamma(xm
+ 1.0);
677 k
= floor(sqrt2xm
* y
+ xm
);
679 t
= 0.9 * (1.0 + y
* y
) * exp(k
* logxm
- gamma(k
+ 1.0) - g
);
685 /*****************************************************************************/
686 /*****************************************************************************/
688 /** Binomial Distribution **/
690 /*****************************************************************************/
691 /*****************************************************************************/
693 static double binomial_cdf(k
, n
, p
)
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;
707 ia
= floor(da
); ib
= floor(db
);
708 betabase(&p
, &da
, &db
, &ia
, &ib
, &dp
);
714 double binomialcdf(k
, n
, p
)
719 return(binomial_cdf((int) floor(k
), n
, p
));
723 int binomialquant(x
, n
, p
)
727 int k
, k1
, k2
, del
, ia
;
728 double m
, s
, p1
, p2
, pk
;
731 checkprob(x
, FALSE
, FALSE
);
734 s
= sqrt(n
* p
* (1 - p
));
735 del
= max(1, (int) (0.2 * s
));
737 if (x
== 0.0) k
= 0.0;
738 else if (x
== 1.0) k
= n
;
739 else k
= m
+ s
* ppnd(x
, &ia
);
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
);
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
);
754 while (k2
- k1
> 1) {
756 pk
= binomial_cdf(k
, n
, p
);
757 if (pk
< x
) { k1
= k
; p1
= pk
; }
758 else { k2
= k
; p2
= pk
; }
763 double binomialpmf(k
, n
, p
)
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
));
779 /* binomial random generator from Numerical Recipes */
780 int binomialrand(n
, pp
)
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
;
789 checkbinomial(n
, pp
);
791 p
= (pp
<= 0.5) ? pp
: 1.0 - pp
;
795 else if (p
== 1.0) k
= n
;
798 for (j
= 0; j
< n
; j
++) if (uni() < p
) k
++;
813 oldg
= gamma(en
+ 1.0);
822 sq
= sqrt(2.0 * am
* pc
);
827 } while (em
< 0.0 || em
>= en
+ 1.0);
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
);
835 if (p
!= pp
) k
= n
- k
;