new randist/binomial.c
James Theiler
jt@lanl.gov
Fri May 9 16:12:00 GMT 2003
Hello,
Awhile back I wrote about how much slower the GSL binomial random
variates were than competing algorithms -- thanks again to John
Pearson for first pointing this out to me.
I have attached new code for binomial variates which incorporates a
state of the art algorithm (called BTPE [*], the same basic algorithm
that is implemented in ranlib, and I am grateful to Dr.
Kachitvichyanukul and Dr. Schmeiser for permission to modify that code
and include it in GSL; also for useful comments and discussion about
the code and random variate generation in general). The new code is
about an order of magnitude faster than our original algorithm, and
several times faster than the Numerical Recipes algorithm. Although
it is essentially the same algorithm that appears in ranlib, I've
managed through various tweaks to produce something that is a little
bit faster than the ranlib code (at least on my machine).
More details are available in the introductory comments to the
source code, but a few points of general interest:
1. Following GSL conventions, this algorithm does not employ static
variables. However, by using static variables to hold intermediate
computations, you can produce significantly (~50%) faster variates in
the situation where you call the function multiple times with the same
parameters. I have to admit that I found this performance gain too
seductive to resist, so I have included in the source code a '#define
USE_XSTATIC' (the XSTATIC isn't supposed to convey my excitement, just
wanted to avoid possible name conflict with USE_STATIC) which is by
default set to zero, but which an individual can set to 1 to get
faster performance.
2. As if that wasn't enough, I also implemented another approach --
one that doesn't use static variables -- for speeding up multiple
calls with the same parameters (for binomial, they are n and p). The
formulation follows that of gsl_ran_discrete (and in fact uses the
gsl_ran_discrete algorithm as well); take a look at the file
binomial_batch.c for details.
3. I have added a new function gsl_ran_binomial_cdf(), which is the
cumulative distribution function (akin to the pdf). I think we should
begin filling out the cdf's for all our nonuniform variates.
4. Just a quick note, I got a substantial performance gain using
gsl_sf_pow_int() in place of pow(), even for very large powers. GSL
developers might keep that in mind when performing x^n type
computations. By the way, I noticed GSL also has a gsl_pow_int()
with essentially the same implementation in the sys/pow_int.c file;
which is preferred?
regards,
jt
[*] Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random
Variate Generation. Communications of the ACM, 31, 2 (February,
1988) 216. There is also an algorithm called PTPE for Poisson
variates; Dr. Schmeiser gave me a fortran implementation and I
have GSL'ified that as well, and will be posting it soon.
---------------------------------------------
James Theiler jt@lanl.gov
MS-B244, NIS-2, LANL tel: 505/665-5682
Los Alamos, NM 87545 fax: 505/665-4414
----- Space and Remote Sensing Sciences -----
-------------- next part --------------
? binomial_batch.c
? discrete-JT.c
? new-binomial.diff
? new-binomial.tar.gz
Index: Makefile.am
===================================================================
RCS file: /cvs/gsl/gsl/randist/Makefile.am,v
retrieving revision 1.44
diff -u -r1.44 Makefile.am
--- Makefile.am 10 Dec 2002 19:06:57 -0000 1.44
+++ Makefile.am 9 May 2003 07:24:20 -0000
@@ -4,7 +4,7 @@
INCLUDES= -I$(top_builddir)
-libgslrandist_la_SOURCES = bernoulli.c beta.c bigauss.c binomial.c cauchy.c chisq.c dirichlet.c discrete.c erlang.c exponential.c exppow.c fdist.c flat.c gamma.c gauss.c gausstail.c geometric.c gumbel.c hyperg.c laplace.c levy.c logarithmic.c logistic.c lognormal.c multinomial.c nbinomial.c pareto.c pascal.c poisson.c rayleigh.c shuffle.c sphere.c tdist.c weibull.c landau.c
+libgslrandist_la_SOURCES = bernoulli.c beta.c bigauss.c binomial.c binomial_batch.c cauchy.c chisq.c dirichlet.c discrete.c erlang.c exponential.c exppow.c fdist.c flat.c gamma.c gauss.c gausstail.c geometric.c gumbel.c hyperg.c laplace.c levy.c logarithmic.c logistic.c lognormal.c multinomial.c nbinomial.c pareto.c pascal.c poisson.c rayleigh.c shuffle.c sphere.c tdist.c weibull.c landau.c
TESTS = test
Index: binomial.c
===================================================================
RCS file: /cvs/gsl/gsl/randist/binomial.c,v
retrieving revision 1.13
diff -u -r1.13 binomial.c
--- binomial.c 24 Apr 2001 17:26:42 -0000 1.13
+++ binomial.c 9 May 2003 07:24:21 -0000
@@ -1,6 +1,6 @@
/* randist/binomial.c
*
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
+ * Copyright (C) 1996-2003 James Theiler, Brian Gough
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -22,52 +22,359 @@
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
+#include <gsl/gsl_sf_pow_int.h>
/* The binomial distribution has the form,
- prob(k) = n!/(k!(n-k)!) * p^k (1-p)^(n-k) for k = 0, 1, ..., n
+ f(x) = n!/(x!(n-x)!) * p^x (1-p)^(n-x) for integer 0 <= x <= n
+ = 0 otherwise
- This is the algorithm from Knuth */
+ This implementation follows the public domain ranlib function
+ "ignbin", the bulk of which is the BTPE (Binomial Triangle
+ Parallelogram Exponential) algorithm introduced in K+S [*], and
+ tries to use more conventional C coding standards. While C
+ ordinarily eschews the goto, and this implementation has reduced
+ the number of goto's substantially, there are still two labels
+ (TryAgain, and Finish) and a number of places where the code uses
+ goto TryAgain, and goto Finish. Since the algorithm involves many
+ accept/reject checks, this is arguably appropriate: the "Finish"
+ corresponds to "accept" and the "TryAgain" to "reject".
+
+ if n is small and/or p is near 0 or near 1 (specifically, if
+ n*min(p,1-p) < SMALL_MEAN), then a different algorithm, called
+ BINV in K+S, is used which has an average runtime that scales
+ linearly with n*min(p,1-p).
+
+ but for larger problems, the BTPE algorithm takes the form of two
+ functions b(x) and t(x) -- "bottom" and "top" -- for which
+ b(x) < f(x)/f(M) < t(x), with M = floor(n*p+p). b(x) defines
+ a triangular region, and t(x) includes a parallelogram and two
+ tails. Details (including a nice drawing) are in K+S.
+
+ [*] Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random
+ Variate Generation. Communications of the ACM, 31, 2 (February,
+ 1988) 216.
+
+ Note, Bruce Schmeiser (personal communication) points out that if
+ you want very fast binomial deviates, and you are happy with
+ approximate results, and/or n and n*p are both large, then you can
+ just use gaussian estimates: mean=n*p, variance=n*p*(1-p).
+
+ Another note, the default form of this version has the same cost
+ per call, whether called once or called many times with the same
+ parameters. By setting the variable USE_XSTATIC to 1, the code
+ will do multiple calls with the same parameters much more cheaply
+ (typically 50% faster).
+
+ Another approach, based on Walker's algorithm as implemented in the
+ gsl_ran_discrete() function is also implemented in binomial_batch.c
+ which provides gsl_ran_binomial_preproc and gsl_ran_binomial_batch.
+ This avoids the problems with static variables; it has a larger
+ setup time but then a faster marginal time once it is set up. You
+ should use this only if you will be running at least a few thousand
+ variates with the same parameters.
+
+ This implementation by James Theiler, April 2003, after obtaining
+ permission -- and some good advice -- from Drs. Kachitvichyanukul
+ and Schmeiser to use their code as a starting point, and then doing
+ a little bit of tweaking.
+
+*/
+
+#define USE_XSTATIC 0 /* Using static variables speeds things up
+ * when there are multiple calls with the
+ * same parameters; but GSL generally eschews
+ * static variables -- not re-entrant, etc.
+ */
+#if USE_XSTATIC
+#define XSTATIC static
+#else
+#define XSTATIC
+#endif
+
+
+#define SMALL_MEAN 14 /* n*p < SMALL_MEAN -> use BINV algorithm */
+ /* ranlib implementation used cutoff=30;
+ * but on my computer 14 works better */
+#define SMALL_MEAN_LARGE_N 110 /* in BINV, do not permit ix too large */
+#define FAR_FROM_MEAN 20 /* ix-n*p larger than this, then "squeeze" */
+ /* ranlib used 20, and this seems to be the
+ * best choice on my machine as well
+ */
unsigned int
-gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n)
+gsl_ran_binomial (const gsl_rng * rng, double pp, unsigned int n)
{
- unsigned int i, a, b, k = 0;
+ int ix; /* return value */
- while (n > 10) /* This parameter is tunable */
- {
- double X;
- a = 1 + (n / 2);
- b = 1 + n - a;
-
- X = gsl_ran_beta (r, (double) a, (double) b);
-
- if (X >= p)
- {
- n = a - 1;
- p /= X;
- }
- else
- {
- k += a;
- n = b - 1;
- p = (p - X) / (1 - X);
- }
+ /* setup */
+ double p = (pp <= 0.5) ? pp : 1-pp; /* p=min(pp,1-pp) */
+ double xnp = n*p;
+ XSTATIC double pprev=-1;
+ XSTATIC unsigned int nprev=0;
+
+ if (n==0) return 0;
+
+ if (xnp < SMALL_MEAN) {
+ /* Inverse cdf logic for small mean
+ * Called BINV in K+S
+ */
+ XSTATIC double fo,s;
+ if (!USE_XSTATIC || ! p==pprev || ! n==nprev) {
+ double q = 1-p;
+ fo = gsl_sf_pow_int(q,n); /* f(x), starting with x=0 */
+ s = p/q;
+ }
+ while (1) {
+ /* This while(1) loop will almost certainly only loop once; but
+ * if u=1 to within a few epsilons of machine precision, then it
+ * is possible for roundoff to prevent the main loop over ix to
+ * achieve its proper value. following the ranlib implementation,
+ * we introduce a check for that situation, and when it occurs,
+ * we just try again.
+ */
+ double f=fo;
+ double u = gsl_rng_uniform(rng);
+ for (ix=0; ix<=SMALL_MEAN_LARGE_N; ++ix) {
+ if(u < f) goto Finish;
+ u -= f;
+ /* Use recursion f(x+1) = f(x)*[(n-x)/(x+1)]*[p/(1-p)] */
+ f *= s*(n-ix)/(ix+1);
+ }
+ /* It should be the case that the 'goto Finish' was encountered
+ * before this point was ever reached. But if we have reached
+ * this point, then roundoff has prevented u from decreasing
+ * all the way to zero. This can happen only if the initial u
+ * was very nearly equal to 1, which is a rare situation. In
+ * that rare situation, we just try again.
+ *
+ * Note, following the ranlib implementation, we loop ix only to
+ * a hardcoded value of SMALL_MEAN_LARGE_N=110; we could have
+ * looped to n, and 99.99...% of the time it won't matter. This
+ * choice, I think is a little more robust against the rare
+ * roundoff error. If n>LARGE_N, then it is technically
+ * possible for ix>LARGE_N, but it is astronomically rare, and
+ * if ix is that large, it is more likely due to roundoff than
+ * probability, so better to nip it at LARGE_N than to take a
+ * chance that roundoff will somehow conspire to produce an even
+ * larger (and more improbable) ix. If n<LARGE_N, then once
+ * ix=n, f=0, and the loop will continue until ix=LARGE_N.
+ */
+ }
+ }
+ else {
+ /* for larger n >= SMALL_MEAN, we invoke the BTPE algorithm */
+
+ int k;
+ double u,v; /* random variates */
+ XSTATIC int m; /* argmax_x f(x) */
+ XSTATIC double p1,p2,p3,p4; /* cumulative area of tri, para, exp tails */
+ XSTATIC double c,xl,xr,lambda_l,lambda_r; /* parameters of para, exp tails */
+ XSTATIC double a,ffm,fm,q,x,xm,xnpq; /* various intermediate values... */
+ if (!USE_XSTATIC || ! p==pprev || ! n==nprev) {
+ q = 1-p;
+ ffm = xnp+p; /* ffm = n*p+p */
+ m = (int)ffm; /* m = int floor[n*p+p] */
+ fm = m; /* fm = double m; */
+ xm = fm+0.5; /* xm = half integer mean (tip of triangle) */
+ xnpq = xnp*q; /* xnpq = n*p*q */
+
+ /* p1: radius of triangle region; since height=1, also: area of region */
+ /* p2: p1 + area of parallelogram region */
+ /* p3: p2 + area of left tail */
+ /* p4: p3 + area of right tail */
+ /* pi/p4: probability of i'th area (i=1,2,3,4) */
+
+ /* Note: magic numbers 2.195, 4.6, 0.134, 20.5, 15.3 */
+ /* These magic numbers are not adjustable...at least not easily! */
+
+ p1 = floor(2.195*sqrt(xnpq)-4.6*q)+0.5;
+
+ /* xl, xr: left and right edges of triangle */
+ xl = xm-p1;
+ xr = xm+p1;
+
+ /* Left tail: t(x) = c*exp(-lambda_l*[xl - (x+0.5)]) */
+ /* Right tail: t(x) = c*exp(-lambda_r*[(x+0.5) - xr]) */
+ c = 0.134+20.5/(15.3+fm);
+ a = (ffm-xl)/(ffm-xl*p); lambda_l = a*(1.0+0.5*a);
+ a = (xr-ffm)/(xr*q); lambda_r = a*(1.0+0.5*a);
+ p2 = p1*(1.0+c+c);
+ p3 = p2+c/lambda_l;
+ p4 = p3+c/lambda_r;
}
- for (i = 0; i < n; i++)
- {
- double u = gsl_rng_uniform (r);
- if (u < p)
- k++;
+ TryAgain:
+ /* generate random variates */
+ u = gsl_rng_uniform(rng)*p4; /* specifies which region: Tri, Par, Tail */
+ v = gsl_rng_uniform(rng);
+
+ if (u > p1) {
+ if (u > p2) {
+ if (u > p3) {
+ /* Right tail */
+ ix = (int)(xr-log(v)/lambda_r);
+ if(ix > n) goto TryAgain;
+ v *= ((u-p3)*lambda_r);
+ }
+ else {
+ /* Left tail */
+ ix = (int)(xl+log(v)/lambda_l);
+ if(ix < 0) goto TryAgain;
+ v *= ((u-p2)*lambda_l);
+ }
+ }
+ else {
+ /* Parallelogram region */
+ x = xl+(u-p1)/c;
+ v = v*c + 1.0 - ((x < xm) ? xm-x : x-xm)/p1;
+ if(v > 1.0 || v <= 0.0) goto TryAgain;
+ ix = (int)x;
+ }
+ }
+ else {
+ /* Triangular region */
+ ix = (int)(xm-p1*v+u);
+ goto Finish;
}
- return k;
+ /* At this point, the goal is to test whether v <= f(x)/f(M)
+ */
+
+ if (0) {
+ /* Here is a direct test. It is a little slower
+ * than the various "squeezing" computations below, but
+ * if things are working, it should give exactly the same
+ * answer (given the same random number seed).
+ */
+ /* Use
+ *
+ * f(x) m!(n-m)!
+ * v <= ---- = -------- * (p/q)^{x-m}
+ * f(M) x!(n-x)!
+ *
+ * ie,
+ *
+ * log(v) <= log(m!) + log((n-m)!) - log(x!) - log((n-x)!) + (x-m)*log(p/q)
+ *
+ */
+ /* Actually, if n<170, then this is sometimes (not always) faster
+ * than the squeezing below, since gsl_sf_lnfact() is based on
+ * precomputed values.
+ */
+ XSTATIC double lnmnm,lnpq;
+ if (!USE_XSTATIC || ! p==pprev || ! n==nprev) {
+ lnmnm = gsl_sf_lnfact(m) + gsl_sf_lnfact(n-m);
+ lnpq = log(p/q);
+ }
+ if (log(v) <= lnmnm -
+ gsl_sf_lnfact(ix) - gsl_sf_lnfact(n-ix) + (ix-m)*lnpq) {
+ goto Finish;
+ } else {
+ goto TryAgain;
+ }
+ }
+
+ /* More efficient determination of whether v < f(x)/f(M) */
+ k = ix > m ? ix-m : m-ix;
+ if (k > FAR_FROM_MEAN) {
+ /* If ix is far from the mean m: k=ABS(ix-m) large
+ */
+ double x1,x2,w1,w2,f1,f2,z1,z2;
+ double alv = log(v);
+ if (k < xnpq/2 - 1) {
+ /* "Squeeze" using upper and lower bounds on log(f(x))
+ * The squeeze condition was derived under the condition k < xnpq/2-1
+ */
+ double amaxp = k/xnpq*((k*(k/3.0+0.625)+(1.0/6.0))/xnpq+0.5);
+ double ynorm = -(k*k/(2.0*xnpq));
+ if(alv < ynorm-amaxp) goto Finish;
+ if(alv > ynorm+amaxp) goto TryAgain;
+ }
+
+#if 1
+ /* Now, again: do the test log(v) vs. log f(x)/f(M) */
+
+ /* The "#define Stirling" below corresponds to the first five
+ * terms in asymptoic formula for
+ * log Gamma (y) - (y-0.5)log(y) + y - 0.5 log(2*pi);
+ * See Abramowitz and Stegun, eq 6.1.40
+ * Here y2 is y1*y1
+ */
+
+#define Stirling(y1,y2) (13860.0-(462.0-(132.0-(99.0-140.0/y2)/y2)/y2)/y2)/y1/166320.0
+
+ /* Note below: two Stirling's are added, and two are subtracted.
+ * In both K+S, and in the ranlib implementation, all four are
+ * added. I (jt) believe that is a mistake -- this has been confirmed
+ * by personal correspondence w/ Dr. Kachitvichyanukul. Note, however,
+ * the corrections are so small, that I couldn't find an example where
+ * it made a difference that could be observed, let alone tested. In
+ * fact, define'ing Stirling to be zero gave identical results!! In
+ * practice, alv is O(1), ranging 0 to -10 or so, while the Stirling
+ * correction is typically O(10^{-5}) ...setting the correction to zero
+ * gives about a 2% performance boost; might as well keep it just to
+ * be pendantic.
+ */
+
+ x1 = ix+1.0; x2 = x1*x1;
+ w1 = n-ix+1.0; w2 = w1*w1;
+ f1 = fm+1.0; f2 = f1*f1;
+ z1 = n+1.0-fm; z2 = z1*z1;
+
+ if(alv <= xm*log(f1/x1)+(n-m+0.5)*log(z1/w1)+(ix-m)*log(w1*p/(x1*q))+
+ Stirling(f1,f2) + Stirling(z1,z2) - Stirling(x1,x2) - Stirling(w1,w2)) {
+ goto Finish;
+ } else {
+ goto TryAgain;
+ }
+#else
+ /* This is equivalent to the above, but is a little (~20%) slower */
+ /* There are five log's vs three above, maybe that's it? */
+ if (alv <= gsl_sf_lnfact(m) + gsl_sf_lnfact(n-m) -
+ gsl_sf_lnfact(ix) - gsl_sf_lnfact(n-ix) + (ix-m)*log(p/q)) {
+ goto Finish;
+ } else {
+ goto TryAgain;
+ }
+#endif
+ }
+ else {
+ /*
+ * If ix near m (ie, |ix-m|<FAR_FROM_MEAN), then do
+ * explicit evaluation using recursion relation for f(x)
+ */
+ double s = p/q;
+ double g = (n+1)*s;
+ double f = 1.0;
+ if(m < ix) {
+ int i;
+ for(i=m+1; i<=ix; i++) {
+ f *= (g/i-s);
+ }
+ }
+ else if(m > ix) {
+ int i;
+ for(i=ix+1; i<=m; i++) {
+ f /= (g/i-s);
+ }
+ }
+ if(v <= f) goto Finish;
+ else goto TryAgain;
+ }
+ }
+
+ Finish:
+#if USE_XSTATIC
+ nprev = n;
+ pprev = p;
+#endif
+ return (pp > 0.5 ) ? n - ix : ix;
}
double
gsl_ran_binomial_pdf (const unsigned int k, const double p,
- const unsigned int n)
+ const unsigned int n)
{
if (k > n)
{
@@ -75,13 +382,64 @@
}
else
{
- double a = k;
- double b = n - k;
double P;
- double Cnk = gsl_sf_choose (n, k) ;
-
- P = Cnk * pow (p, a) * pow (1 - p, b);
+ if (n>170)
+ {
+ /* Hardcoded 170 corresponds to the largest n!
+ * that can be represented as a double precision
+ * number, before overflow. In that case, use
+ * logarithms.
+ */
+ double ln_Cnk = gsl_sf_lnchoose (n, k);
+ P = ln_Cnk + k*log(p) + (n-k)*log(1-p);
+ P = exp(P);
+ }
+ else
+ {
+ double Cnk = gsl_sf_choose (n, k);
+ P = Cnk * gsl_sf_pow_int(p,k) * gsl_sf_pow_int(1-p,n-k);
+ }
return P;
}
+}
+double
+gsl_ran_binomial_cdf (const unsigned int k, const double p,
+ const unsigned int n)
+{
+ /* cdf returns pdf(0)+pdf(1)+...+pdf(k) */
+ /* The cdf is computed by first computing pdf(k) and then
+ * using a recursion scheme to compute values of pdf(k-1),
+ * pdf(k-2), ... , pdf(0), summing them up as it goes along,
+ * and bailing out early if the pdf reaches numerical zero.
+ *
+ * Note that if k>p*n, then we use a slightly different strategy
+ * we use the recursion to compute pdf(k+1), pdf(k+2), ..., until
+ * pdf(n), summing those as it goes along. Then this is 1-cdf.
+ */
+ double Ck=0;
+ if (p>0) /* if p==0, then return Ck=0; */
+ {
+ double s = p/(1-p);
+ double Pk;
+ int i;
+ if (k < n*p)
+ {
+ Ck = Pk = gsl_ran_binomial_pdf (k,p,n);
+ for (i=k; Pk>0; --i) {
+ Pk = Pk * i/(s*(n-i+1));
+ Ck += Pk;
+ }
+ }
+ else
+ {
+ Ck = Pk = gsl_ran_binomial_pdf(k+1,p,n);
+ for (i=k+1; Pk>0; ++i) {
+ Pk = Pk * s*(n-i)/(i+1);
+ Ck += Pk;
+ }
+ Ck = 1-Ck;
+ }
+ }
+ return Ck;
}
Index: gsl_randist.h
===================================================================
RCS file: /cvs/gsl/gsl/randist/gsl_randist.h,v
retrieving revision 1.39
diff -u -r1.39 gsl_randist.h
--- gsl_randist.h 10 Dec 2002 19:06:57 -0000 1.39
+++ gsl_randist.h 9 May 2003 07:24:21 -0000
@@ -41,6 +41,7 @@
unsigned int gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n);
double gsl_ran_binomial_pdf (const unsigned int k, const double p, const unsigned int n);
+double gsl_ran_binomial_cdf (const unsigned int k, const double p, const unsigned int n);
double gsl_ran_exponential (const gsl_rng * r, const double mu);
double gsl_ran_exponential_pdf (const double x, const double mu);
@@ -173,6 +174,17 @@
void gsl_ran_discrete_free(gsl_ran_discrete_t *g);
size_t gsl_ran_discrete (const gsl_rng *r, const gsl_ran_discrete_t *g);
double gsl_ran_discrete_pdf (size_t k, const gsl_ran_discrete_t *g);
+
+typedef struct {
+ unsigned int n;
+ double p;
+ unsigned int noff; /* offset, k<noff => P[k]=0 */
+ gsl_ran_discrete_t *discrete;
+} gsl_ran_binomial_t;
+
+void gsl_ran_binomial_free(gsl_ran_binomial_t *binom);
+unsigned int gsl_ran_binomial_batch(gsl_rng *rng, gsl_ran_binomial_t *binom);
+gsl_ran_binomial_t *gsl_ran_binomial_preproc(double p, unsigned int n);
__END_DECLS
Index: test.c
===================================================================
RCS file: /cvs/gsl/gsl/randist/test.c,v
retrieving revision 1.38
diff -u -r1.38 test.c
--- test.c 9 Feb 2003 23:23:08 -0000 1.38
+++ test.c 9 May 2003 07:24:21 -0000
@@ -29,7 +29,7 @@
#define N 100000
-/* Convient test dimension for multivariant distributions */
+/* Convient test dimension for multivariate distributions */
#define MULTI_DIM 10
@@ -38,7 +38,8 @@
void testPDF (double (*f) (void), double (*pdf) (double), const char *name);
void testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
const char *name);
-
+void testDiscreteCDF (double (*pdf) (unsigned int),
+ double (*cdf) (unsigned int), const char *name);
void test_shuffle (void);
void test_choose (void);
double test_beta (void);
@@ -47,8 +48,16 @@
double test_bernoulli_pdf (unsigned int n);
double test_binomial (void);
double test_binomial_pdf (unsigned int n);
+double test_binomial_cdf (unsigned int n);
double test_binomial_large (void);
double test_binomial_large_pdf (unsigned int n);
+double test_binomial_large_cdf (unsigned int n);
+double test_binomial_huge (void);
+double test_binomial_huge_pdf (unsigned int n);
+double test_binomial_huge_cdf (unsigned int n);
+double test_binomial_batch (void);
+double test_binomial_batch_pdf (unsigned int n);
+double test_binomial_batch_cdf (unsigned int n);
double test_cauchy (void);
double test_cauchy_pdf (double x);
double test_chisq (void);
@@ -212,6 +221,7 @@
#define FUNC(x) test_ ## x, "test gsl_ran_" #x
#define FUNC2(x) test_ ## x, test_ ## x ## _pdf, "test gsl_ran_" #x
+#define FUNC2CDF(x) test_ ## x ## _pdf, test_ ## x ## _cdf, "test PDF/CDF " #x
test_shuffle ();
test_choose ();
@@ -308,6 +318,8 @@
testDiscretePDF (FUNC2 (bernoulli));
testDiscretePDF (FUNC2 (binomial));
testDiscretePDF (FUNC2 (binomial_large));
+ testDiscretePDF (FUNC2 (binomial_huge));
+ testDiscretePDF (FUNC2 (binomial_batch));
testDiscretePDF (FUNC2 (geometric));
testDiscretePDF (FUNC2 (geometric1));
testDiscretePDF (FUNC2 (hypergeometric1));
@@ -322,6 +334,10 @@
testDiscretePDF (FUNC2 (negative_binomial));
testDiscretePDF (FUNC2 (pascal));
+ testDiscreteCDF (FUNC2CDF (binomial));
+ testDiscreteCDF (FUNC2CDF (binomial_large));
+ testDiscreteCDF (FUNC2CDF (binomial_huge));
+
exit (gsl_test_summary ());
}
@@ -550,6 +566,37 @@
name, 0, BINS);
}
+void
+testDiscreteCDF (double (*pdf) (unsigned int),
+ double (*cdf) (unsigned int), const char *name)
+{
+ double c[BINS], p[BINS];
+ unsigned int i;
+ int status = 0, status_i = 0;
+ double count=0;
+
+
+ for (i = 0; i < BINS; i++)
+ {
+ double dc=0;
+ p[i] = pdf (i);
+ c[i] = cdf (i);
+ count += p[i];
+ dc = count - c[i];
+ status_i = (dc*dc > 1.0e-10);
+ status |= status_i;
+ if (status_i)
+ {
+ gsl_test (status_i, "%s i=%d (%g observed vs %g expected)",
+ name, i, count , c[i]);
+ break;
+ }
+ }
+
+ if (status == 0)
+ gsl_test (status, "%s, comparing cdf against pdf over range [%d,%d) ",
+ name, 0, BINS);
+}
double
@@ -588,6 +635,11 @@
{
return gsl_ran_binomial_pdf (n, 0.3, 5);
}
+double
+test_binomial_cdf (unsigned int n)
+{
+ return gsl_ran_binomial_cdf (n, 0.3, 5);
+}
double
test_binomial_large (void)
@@ -599,6 +651,47 @@
test_binomial_large_pdf (unsigned int n)
{
return gsl_ran_binomial_pdf (n, 0.3, 55);
+}
+double
+test_binomial_large_cdf (unsigned int n)
+{
+ return gsl_ran_binomial_cdf (n, 0.3, 55);
+}
+
+double
+test_binomial_huge (void)
+{
+ return gsl_ran_binomial (r_global, 0.3, 5500);
+}
+
+double
+test_binomial_huge_pdf (unsigned int n)
+{
+ return gsl_ran_binomial_pdf (n, 0.3, 5500);
+}
+double
+test_binomial_huge_cdf (unsigned int n)
+{
+ return gsl_ran_binomial_cdf (n, 0.3, 5500);
+}
+
+static gsl_ran_binomial_t *binom=NULL;
+double
+test_binomial_batch (void)
+{
+ if (binom == NULL) binom = gsl_ran_binomial_preproc(0.3,550);
+ return gsl_ran_binomial_batch (r_global, binom);
+}
+
+double
+test_binomial_batch_pdf (unsigned int n)
+{
+ return gsl_ran_binomial_pdf (n, 0.3, 550);
+}
+double
+test_binomial_batch_cdf (unsigned int n)
+{
+ return gsl_ran_binomial_cdf (n, 0.3, 550);
}
double
-------------- next part --------------
A non-text attachment was scrubbed...
Name: new-binomial.tar.gz
Type: application/x-gunzip
Size: 13882 bytes
Desc:
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20030509/fd1924f8/attachment.bin>
More information about the Gsl-discuss
mailing list