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(),
15 /* forward declaration */
16 extern double tdens();
19 #define PI 3.14159265358979323846
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
40 return(sin(x
) / cos(x
));
43 double floor(double x
)
47 return((dx
<= x
) ? dx
: dx
- 1.0);
55 /* do nothing for now */
62 xlfail("non-positive gamma or beta exponent");
70 xlfail("non-positive degrees of freedom");
75 checkprob(double p
, long zerostrict
, long onestrict
)
78 if (p
<= 0.0) xlfail("non-positive probability argument");
80 if (p
< 0.0) xlfail("negative probability argument");
83 if (p
>= 1.0) xlfail("probability argument not less than one");
85 if (p
> 1.0) xlfail("probability argument greater than one");
92 if (r
< -1 || r
> 1) {
93 xlfail("correlation out of range");
98 checkpoisson(double L
)
101 xlfail("negative Poisson mean");
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 /*****************************************************************************/
115 /** Uniform Distribution **/
117 /*****************************************************************************/
118 /*****************************************************************************/
120 /* uniform generator - avoids zero and one */
126 } while ((u
<= 0.0) || (u
>= 1.0));
130 /*****************************************************************************/
131 /*****************************************************************************/
133 /** Normal Distribution **/
135 /*****************************************************************************/
136 /*****************************************************************************/
138 /* standard normal cdf */
139 double normalcdf(double x
)
146 /* standard normal quantile function */
147 double normalquant(double p
)
152 checkprob(p
, TRUE
, TRUE
);
158 /* standard normal density */
159 double normaldens(double x
)
161 return(exp(- 0.5 * x
* x
) / sqrt(2.0 * PI
));
164 /* standard normal generator */
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 */
176 v
= c
* (2 * u1
- 1);
179 } while(y
> (1 - u
) && y
> - log(u
));
183 /* bivariate normal cdf */
184 double bnormcdf(double x
, double y
, double r
)
187 return(bivnor(-x
, -y
, r
));
190 /*****************************************************************************/
191 /*****************************************************************************/
193 /** Cauchy Distribution **/
195 /*****************************************************************************/
196 /*****************************************************************************/
201 return((atan(dx
) + PI
/ 2) / PI
);
205 cauchyquant(double p
)
207 checkprob(p
, TRUE
, TRUE
);
208 return(tan(PI
* (p
- 0.5)));
212 cauchydens(double dx
)
214 return(tdens(dx
, 1.0));
217 /* cauchy generator */
221 double u1
, u2
, v1
, v2
;
223 /* ratio of uniforms on half disk */
229 } while(v1
* v1
+ v2
* v2
> 1.0);
233 /*****************************************************************************/
235 /** Gamma Distribution **/
237 /*****************************************************************************/
240 gammacdf(double x
, double a
)
245 if (x
<= 0.0) p
= 0.0;
246 else gammabase(&x
, &a
, &p
);
251 gammaquant(double p
, double a
)
257 checkprob(p
, FALSE
, TRUE
);
258 x
= ppgamma(p
, a
, &flag
);
264 double gammadens(double x
, double a
)
272 dens
= exp(log(x
) * (a
- 1) - x
- gamma(a
));
277 /* gamma generator */
281 double x
, u0
, u1
, u2
, v
, w
, c
, c1
, c2
, c3
, c4
, c5
;
282 static double e
= -1.0;
286 if (e
< 0.0) e
= exp(1.0);
289 /* Ahrens and Dieter algorithm */
298 if (u1
<= exp(-x
)) done
= TRUE
;
301 x
= -log((c
- v
) / a
);
302 if (x
> 0.0 && u1
< exp((a
- 1.0) * log(x
))) done
= TRUE
;
306 else if (a
== 1.0) x
= -log(unirand());
308 /* Cheng and Feast algorithm */
310 c2
= (a
- 1.0 / (6.0 * a
)) / c1
;
312 c4
= 2.0 / (a
- 1.0) + 2.0;
318 if (a
> 2.5) u1
= u2
+ c5
* (1.0 - 1.86 * u1
);
319 } while (u1
<= 0.0 || u1
>= 1.0);
321 } while ((c3
* u1
+ w
+ 1.0/w
) > c4
&& (c3
* log(u1
) - log(w
) + w
) > 1.0);
327 /*****************************************************************************/
329 /** Chi-Square Distribution **/
331 /*****************************************************************************/
333 double chisqcdf(double x
, double df
)
343 gammabase(&x
, &a
, &p
);
349 chisqquant(double p
, double df
)
355 checkprob(p
, FALSE
, TRUE
);
357 x
= 2.0 * ppgamma(p
, a
, &flag
);
363 chisqdens(double dx
, double da
)
368 return(0.5 * gammadens(dx
, da
));
375 return(2.0 * gammarand(df
/ 2.0));
378 /*****************************************************************************/
380 /** Beta Distribution **/
382 /*****************************************************************************/
385 betacdf(double x
, double a
, double b
)
390 checkexp(a
); checkexp(b
);
398 betabase(&x
, &a
, &b
, &ia
, &ib
, &p
);
405 betaquant(double p
, double a
, double b
)
412 checkprob(p
, FALSE
, FALSE
);
413 x
= ppbeta(p
, a
, b
, &flag
);
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 */
425 lbeta
= gamma(da
) + gamma(db
) - gamma(da
+ db
);
432 betadens(double x
, double a
, double b
)
438 if (x
<= 0.0 || x
>= 1.0) {
441 dens
= exp(log(x
) * (a
- 1) + log(1 - x
) * (b
- 1) - logbeta(a
, b
));
448 betarand(double a
, double b
)
459 /*****************************************************************************/
461 /** t Distribution **/
463 /*****************************************************************************/
467 tcdf(double x
, double df
)
472 studentbase(&x
, &df
, &p
);
477 tquant(double p
, double df
)
483 checkprob(p
, TRUE
, TRUE
);
484 x
= ppstudent(p
, df
, &flag
);
490 tdens(double x
, double 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
));
501 double trand(double df
)
504 return(normalrand() / sqrt(chisqrand(df
) / df
));
507 /*****************************************************************************/
509 /** F Distribution **/
511 /*****************************************************************************/
513 double fcdf(double x
, double ndf
, double ddf
)
517 checkdf(ndf
); checkdf(ddf
);
524 p
= 1.0 - betacdf(x
, a
, b
);
529 double fquant(double p
, double ndf
, double ddf
)
534 checkdf(ndf
); checkdf(ddf
);
535 checkprob(p
, FALSE
, TRUE
);
542 x
= ppbeta(p
, a
, b
, &flag
);
544 x
= a
* (1.0 / x
- 1.0) / b
;
550 fdens(double dx
, double da
, double db
)
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
));
568 double frand(double ndf
, double ddf
)
572 return((ddf
* chisqrand(ndf
)) / (ndf
* chisqrand(ddf
)));
575 /*****************************************************************************/
577 /** Poisson Distribution **/
579 /*****************************************************************************/
582 poisson_cdf(long k
, double L
)
590 dp
= (k
< 0) ? 0.0 : 1.0;
593 gammabase(&L
, &dx
, &dp
);
601 poissoncdf(double k
, double L
)
604 return(poisson_cdf((long) floor(k
), L
));
608 poissonquant(double x
, double L
)
610 long k
, k1
, k2
, del
, ia
;
611 double m
, s
, p1
, p2
, pk
;
614 checkprob(x
, FALSE
, TRUE
);
617 del
= max(1, (long) (0.2 * s
));
622 k
= m
+ s
* ppnd(x
, &ia
);
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
) {
638 p2
= poisson_cdf(k2
, L
);
641 while (k2
- k1
> 1) {
643 pk
= poisson_cdf(k
, L
);
654 poissonpmf(long k
, double L
)
661 dp
= (k
== 0) ? 1.0 : 0.0;
666 dp
= exp(dx
* log(L
) - L
- gamma(dx
+ 1.0));
672 /* poisson random generator from Numerical Recipes */
673 long poissonrand(double xm
)
675 static double sqrt2xm
, logxm
, expxm
, g
, oldxm
= -1.0;
681 if (xm
!= oldxm
) { expxm
= exp(-xm
); oldxm
= xm
; }
691 sqrt2xm
= sqrt(2.0 * xm
);
693 g
= xm
* logxm
- gamma(xm
+ 1.0);
698 k
= floor(sqrt2xm
* y
+ xm
);
700 t
= 0.9 * (1.0 + y
* y
) * exp(k
* logxm
- gamma(k
+ 1.0) - g
);
706 /*****************************************************************************/
708 /** Binomial Distribution **/
710 /*****************************************************************************/
713 binomial_cdf(long k
, long n
, double p
)
725 dp
= (k
< 0) ? 0.0 : 1.0;
728 dp
= (k
< n
) ? 0.0 : 1.0;
734 betabase(&p
, &da
, &db
, &ia
, &ib
, &dp
);
744 binomialcdf(double k
, long n
, double p
)
747 return(binomial_cdf((long) floor(k
), n
, p
));
751 binomialquant(double x
, long n
, double p
)
753 long k
, k1
, k2
, del
, ia
;
754 double m
, s
, p1
, p2
, pk
;
757 checkprob(x
, FALSE
, FALSE
);
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
);
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
);
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) {
782 pk
= binomial_cdf(k
, n
, p
);
783 if (pk
< x
) { k1
= k
; p1
= pk
; }
784 else { k2
= k
; p2
= pk
; }
789 double binomialpmf(long k
, long n
, double p
)
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
));
803 /* binomial random generator from Numerical Recipes */
804 long binomialrand(long n
, double pp
)
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
;
817 else if (p
== 1.0) k
= n
;
820 for (j
= 0; j
< n
; j
++) if (uni() < p
) k
++;
835 oldg
= gamma(en
+ 1.0);
844 sq
= sqrt(2.0 * am
* pc
);
849 } while (em
< 0.0 || em
>= en
+ 1.0);
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
);
857 if (p
!= pp
) k
= n
- k
;