CDF diff

Jason Hooper Stover jason@sakla.net
Wed Jan 15 03:31:00 GMT 2003


Martin,

Thanks for the code, I'll work on adding it to the cdf directory.

> 
> A few related questions and suggestions: could somebody explain the
> rationale for the implementation of the Erlang distribution?  I
> thought it was a discrete distribution, but maybe I'm confused.

According to Kendall&Stuart(&Ord), the Erlang distribution is
another name for the Gamma distribution. They say "Erlang distribution"
is the name sometimes used in the queueing theory literature.

> 
> Second, it might be useful to have a functions that compute a definite
> integral over a pdf with lower and upper bounds supplied by arguments.
> This could be in addition to or in place of the _cdf functions.  The
> reason for this is (discrete) distributions with no closed-form
> formula for the cumulative density/mass.  So the design questions is,
> should one have functions of the form foo_cdf(x,...) where x is the
> upper bound of the sum/integral, or rather foo_cdf(a,b,...) that
> return \sum_{x=a}^b foo_pdf(x,...) (resp. \int)?

The cdf module computes functions of the form Pr(Z<x) and
Pr(Z>x) for a random variable Z, so the argument to the functions
is either the upper or lower bound of the integral. 
I hadn't planned to compute Pr(a<Z<b) directly since such
a function isn't a "cumulative distribution function" 
according to the common definition. Also, users can 
compute Pr(a<Z<b) as Pr(Z<b)-Pr(Z<a).

> 
> Finally, would it be worth the added overhead to have separate types
> for the parameters of each distribution?  For example, the current
> gamma_pdf takes a shape parameter a and a scale parameter b; one could
> put them in a struct and have several constructors (one in terms of
> shape and scale, one in terms of shape in inverse scale, one in terms
> of sample moments, etc.) that return such a struct, which could then
> be passed to the _pdf(), _cdf() and variate functions.  I'd shy away
> from using macros for this, since the macro could only be invoked
> inside the argument list of a function call expression.

My initial reaction is that the added overhead wouldn't be worth the
effort for a library like GSL. E.g., such an approach might make sense 
for a higher-level program that runs statistical
routines for the exponential family. But I could be wrong; maybe it
would make sense in GSL. If you have some coded examples of this 
kind of thing, send them along.

-Jason

> 
> Thanks,
> - martin

> Index: beta.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/beta.c,v
> retrieving revision 1.10
> diff -c -r1.10 beta.c
> *** beta.c	23 Apr 2001 09:38:28 -0000	1.10
> --- beta.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 58,60 ****
> --- 58,77 ----
>         return p;
>       }
>   }
> + 
> + double
> + gsl_ran_beta_cdf (const double x, const double a, const double b)
> + {
> +   if (x <= 0)
> +     {
> +       return 0;
> +     }
> +   else if (x >= 1)
> +     {
> +       return 1;
> +     }
> +   else
> +     {
> +       return gsl_sf_beta_inc (a, b, x);
> +     }
> + }
> Index: binomial.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/binomial.c,v
> retrieving revision 1.13
> diff -c -r1.13 binomial.c
> *** binomial.c	24 Apr 2001 17:26:42 -0000	1.13
> --- binomial.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 85,87 ****
> --- 85,101 ----
>         return P;
>       }
>   }
> + 
> + double
> + gsl_ran_binomial_cdf (const unsigned int k, const double p, 
> + 		      const unsigned int n)
> + {
> +   if (k >= n)
> +     {
> +       return 1;
> +     }
> +   else 
> +     {
> +       return 1 - gsl_sf_beta_inc (k+1, n-k, p);
> +     }
> + }
> Index: cauchy.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/cauchy.c,v
> retrieving revision 1.10
> diff -c -r1.10 cauchy.c
> *** cauchy.c	17 Apr 2001 19:35:56 -0000	1.10
> --- cauchy.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 49,51 ****
> --- 49,57 ----
>     double p = (1 / (M_PI * a)) / (1 + u * u);
>     return p;
>   }
> + 
> + double
> + gsl_ran_cauchy_cdf (const double x, const double a)
> + {
> +   return 0.5 + M_1_PI * atan (x/a);
> + }
> Index: chisq.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/chisq.c,v
> retrieving revision 1.11
> diff -c -r1.11 chisq.c
> *** chisq.c	23 Apr 2001 09:38:28 -0000	1.11
> --- chisq.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 19,24 ****
> --- 19,25 ----
>   
>   #include <config.h>
>   #include <math.h>
> + #include <gsl/gsl_sf_erf.h>
>   #include <gsl/gsl_sf_gamma.h>
>   #include <gsl/gsl_rng.h>
>   #include <gsl/gsl_randist.h>
> ***************
> *** 50,54 ****
> --- 51,74 ----
>         
>         p = exp ((nu / 2 - 1) * log (x/2) - x/2 - lngamma) / 2;
>         return p;
> +     }
> + }
> + 
> + double
> + gsl_ran_chisq_cdf (const double x, const double nu)
> + {
> +   if (x <= 0)
> +     {
> +       return 0;
> +     }
> +   else if (0 && nu == 1)
> +     {
> +       /* It might make sense to treat the case of a single Gaussian
> + 	 specially. */
> +       return gsl_sf_erf ( sqrt (x/2) );
> +     }
> +   else
> +     {
> +       return gsl_sf_gamma_inc_P (nu/2, x/2);
>       }
>   }
> Index: exponential.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/exponential.c,v
> retrieving revision 1.11
> diff -c -r1.11 exponential.c
> *** exponential.c	4 May 2000 11:25:06 -0000	1.11
> --- exponential.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 19,24 ****
> --- 19,25 ----
>   
>   #include <config.h>
>   #include <math.h>
> + #include <gsl/gsl_math.h>
>   #include <gsl/gsl_rng.h>
>   #include <gsl/gsl_randist.h>
>   
> ***************
> *** 45,52 ****
>       }
>     else
>       {
> !       double p = exp (-x/mu)/mu;
>         
>         return p;
>       }
>   }
> --- 46,67 ----
>       }
>     else
>       {
> !       double p = exp (-x/mu) / mu;
>         
>         return p;
> +     }
> + }
> + 
> + double
> + gsl_ran_exponential_cdf (const double x, const double mu)
> + {
> +   if (x <= 0)
> +     {
> +       return 0 ;
> +     }
> +   else
> +     {
> +       /* return 1 - exp (-x/mu); */
> +       return fabs ( gsl_expm1 (-x/mu) );
>       }
>   }
> Index: gamma.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/gamma.c,v
> retrieving revision 1.20
> diff -c -r1.20 gamma.c
> *** gamma.c	23 Apr 2001 09:38:28 -0000	1.20
> --- gamma.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 138,143 ****
> --- 138,146 ----
>     return x;
>   }
>   
> + #define GSL_RAN_GAMMA_SHAPE_SCALE(shape,scale) (shape), 1/(scale)
> + #define GSL_RAN_GAMMA_SHAPE_INVSCALE(shape,invscale) (shape), (invscale)
> + 
>   double
>   gsl_ran_gamma_pdf (const double x, const double a, const double b)
>   {
> ***************
> *** 156,166 ****
>       {
>         return exp(-x/b)/b ;
>       }
> !   else 
>       {
>         double p;
>         double lngamma = gsl_sf_lngamma (a);
>         p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b;
>         return p;
>       }
>   }
> --- 159,187 ----
>       {
>         return exp(-x/b)/b ;
>       }
> !   else
>       {
>         double p;
>         double lngamma = gsl_sf_lngamma (a);
>         p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b;
>         return p;
> +     }
> + }
> + 
> + double
> + gsl_ran_gamma_cdf (const double x, const double a, const double b)
> + {
> +   if (x <= 0)
> +     {
> +       return 0 ;
> +     }
> +   else if (a == 1)
> +     /* exponential_cdf, see exponential.c */
> +     {
> +       return fabs ( gsl_expm1 (-x/b) );
> +     }
> +   else
> +     {
> +       return gsl_sf_gamma_inc_P (a, x/b);
>       }
>   }
> Index: gauss.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/gauss.c,v
> retrieving revision 1.19
> diff -c -r1.19 gauss.c
> *** gauss.c	4 Jun 2002 22:04:43 -0000	1.19
> --- gauss.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 22,27 ****
> --- 22,28 ----
>   #include <gsl/gsl_math.h>
>   #include <gsl/gsl_rng.h>
>   #include <gsl/gsl_randist.h>
> + #include <gsl/gsl_sf_erf.h>
>   
>   /* Of the two methods provided below, I think the Polar method is more
>    * efficient, but only when you are actually producing two random
> ***************
> *** 89,94 ****
> --- 90,102 ----
>     double u = x / fabs (sigma);
>     double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
>     return p;
> + }
> + 
> + double
> + gsl_ran_gaussian_cdf (const double x, const double sigma)
> + {
> +   double u = x / fabs (sigma);
> +   return (1 + gsl_sf_erf (M_SQRT1_2 * u)) / 2;
>   }
>   
>   double
> Index: geometric.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/geometric.c,v
> retrieving revision 1.9
> diff -c -r1.9 geometric.c
> *** geometric.c	4 May 2000 11:25:06 -0000	1.9
> --- geometric.c	10 Jan 2003 18:09:47 -0000
> ***************
> *** 19,24 ****
> --- 19,25 ----
>   
>   #include <config.h>
>   #include <math.h>
> + #include <gsl/gsl_math.h>
>   #include <gsl/gsl_rng.h>
>   #include <gsl/gsl_randist.h>
>   
> ***************
> *** 52,67 ****
>   gsl_ran_geometric_pdf (const unsigned int k, const double p)
>   {
>     if (k == 0)
> !     {
> !       return 0 ;
> !     }
> !   else if (k == 1)
> !     {
> !       return p ;
> !     }
>     else
> !     {
> !       double P = p * pow (1 - p, k - 1.0);
> !       return P;
> !     }
>   }
> --- 53,68 ----
>   gsl_ran_geometric_pdf (const unsigned int k, const double p)
>   {
>     if (k == 0)
> !       return 0;
>     else
> !     return p * gsl_pow_int (1-p, k-1);
> ! }
> ! 
> ! double
> ! gsl_ran_geometric_cdf (const unsigned int k, const double p)
> ! {
> !   if (k == 0)
> !     return 0;
> !   else
> !     return 1 - gsl_pow_int (1-p, k);
>   }
> Index: gsl_randist.h
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/gsl_randist.h,v
> retrieving revision 1.39
> diff -c -r1.39 gsl_randist.h
> *** gsl_randist.h	10 Dec 2002 19:06:57 -0000	1.39
> --- gsl_randist.h	10 Jan 2003 18:09:48 -0000
> ***************
> *** 38,58 ****
> --- 38,63 ----
>   
>   double gsl_ran_beta (const gsl_rng * r, const double a, const double b);
>   double gsl_ran_beta_pdf (const double x, const double a, const double b);
> + double gsl_ran_beta_cdf (const double x, const double a, const double b);
>   
>   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);
> + double gsl_ran_exponential_cdf (const double x, const double mu);
>   
>   double gsl_ran_exppow (const gsl_rng * r, const double a, const double b);
>   double gsl_ran_exppow_pdf (const double x, const double a, const double b);
>   
>   double gsl_ran_cauchy (const gsl_rng * r, const double a);
>   double gsl_ran_cauchy_pdf (const double x, const double a);
> + double gsl_ran_cauchy_cdf (const double x, const double a);
>   
>   double gsl_ran_chisq (const gsl_rng * r, const double nu);
>   double gsl_ran_chisq_pdf (const double x, const double nu);
> + double gsl_ran_chisq_cdf (const double x, const double nu);
>   
>   void gsl_ran_dirichlet (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);
>   double gsl_ran_dirichlet_pdf (const size_t K, const double alpha[], const double theta[]);
> ***************
> *** 70,79 ****
> --- 75,86 ----
>   double gsl_ran_gamma (const gsl_rng * r, const double a, const double b);
>   double gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a);
>   double gsl_ran_gamma_pdf (const double x, const double a, const double b);
> + double gsl_ran_gamma_cdf (const double x, const double a, const double b);
>   
>   double gsl_ran_gaussian (const gsl_rng * r, const double sigma);
>   double gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma);
>   double gsl_ran_gaussian_pdf (const double x, const double sigma);
> + double gsl_ran_gaussian_cdf (const double x, const double sigma);
>   
>   double gsl_ran_ugaussian (const gsl_rng * r);
>   double gsl_ran_ugaussian_ratio_method (const gsl_rng * r);
> ***************
> *** 93,98 ****
> --- 100,106 ----
>   
>   unsigned int gsl_ran_geometric (const gsl_rng * r, const double p);
>   double gsl_ran_geometric_pdf (const unsigned int k, const double p);
> + double gsl_ran_geometric_cdf (const unsigned int k, const double p);
>   
>   unsigned int gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2, unsigned int t);
>   double gsl_ran_hypergeometric_pdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);
> ***************
> *** 105,110 ****
> --- 113,119 ----
>   
>   double gsl_ran_logistic (const gsl_rng * r, const double a);
>   double gsl_ran_logistic_pdf (const double x, const double a);
> + double gsl_ran_logistic_cdf (const double x, const double a);
>   
>   double gsl_ran_lognormal (const gsl_rng * r, const double zeta, const double sigma);
>   double gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma);
> ***************
> *** 129,134 ****
> --- 138,144 ----
>   
>   double gsl_ran_pareto (const gsl_rng * r, double a, const double b);
>   double gsl_ran_pareto_pdf (const double x, const double a, const double b);
> + double gsl_ran_pareto_cdf (const double x, const double a, const double b);
>   
>   unsigned int gsl_ran_poisson (const gsl_rng * r, double mu);
>   void gsl_ran_poisson_array (const gsl_rng * r, size_t n, unsigned int array[],
> ***************
> *** 143,148 ****
> --- 153,159 ----
>   
>   double gsl_ran_tdist (const gsl_rng * r, const double nu);
>   double gsl_ran_tdist_pdf (const double x, const double nu);
> + double gsl_ran_tdist_cdf (const double x, const double nu);
>   
>   double gsl_ran_laplace (const gsl_rng * r, const double a);
>   double gsl_ran_laplace_pdf (const double x, const double a);
> Index: logistic.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/logistic.c,v
> retrieving revision 1.11
> diff -c -r1.11 logistic.c
> *** logistic.c	17 Apr 2001 19:35:56 -0000	1.11
> --- logistic.c	10 Jan 2003 18:09:48 -0000
> ***************
> *** 51,53 ****
> --- 51,59 ----
>     double p = u / (fabs(a) * (1 + u) * (1 + u));
>     return p;
>   }
> + 
> + double
> + gsl_ran_logistic_cdf (const double x, const double a)
> + {
> +   return 1 / (1 + exp(-x/fabs(a)));
> + }
> Index: pareto.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/pareto.c,v
> retrieving revision 1.9
> diff -c -r1.9 pareto.c
> *** pareto.c	21 Sep 2000 19:19:21 -0000	1.9
> --- pareto.c	10 Jan 2003 18:09:48 -0000
> ***************
> *** 52,54 ****
> --- 52,59 ----
>       }
>   }
>   
> + double
> + gsl_ran_pareto_cdf (const double x, const double a, const double b)
> + {
> +   return (x >= b)? 1 - pow (b/x, a) : 0;
> + }
> Index: tdist.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/tdist.c,v
> retrieving revision 1.12
> diff -c -r1.12 tdist.c
> *** tdist.c	23 Apr 2001 09:38:28 -0000	1.12
> --- tdist.c	10 Jan 2003 18:09:48 -0000
> ***************
> *** 77,80 ****
>     return p;
>   }
>   
> ! 
> --- 77,93 ----
>     return p;
>   }
>   
> ! double
> ! gsl_ran_tdist_cdf (const double x, const double nu)
> ! {
> !   if (x == 0)
> !     {
> !       return 0.5;
> !     }
> !   else {
> !     double d = 0.5 * gsl_sf_beta_inc (0.5*nu, 0.5, nu/(nu + x*x));
> !     /* d as a function of x is symmetric about zero, but has a maximum
> !        there; adjust as follows(?): */
> !     return (x < 0)? d : 1 - d;
> !   }
> ! }
> Index: test.c
> ===================================================================
> RCS file: /cvs/gsl/gsl/randist/test.c,v
> retrieving revision 1.37
> diff -c -r1.37 test.c
> *** test.c	10 Dec 2002 19:06:58 -0000	1.37
> --- test.c	10 Jan 2003 18:09:48 -0000
> ***************
> *** 35,58 ****
>   
>   void testMoments (double (*f) (void), const char *name,
>   		  double a, double b, double p);
> ! void testPDF (double (*f) (void), double (*pdf) (double), const char *name);
>   void testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
> ! 		      const char *name);
>   
>   void test_shuffle (void);
>   void test_choose (void);
>   double test_beta (void);
>   double test_beta_pdf (double x);
>   double test_bernoulli (void);
>   double test_bernoulli_pdf (unsigned int n);
>   double test_binomial (void);
>   double test_binomial_pdf (unsigned int n);
>   double test_binomial_large (void);
>   double test_binomial_large_pdf (unsigned int n);
>   double test_cauchy (void);
>   double test_cauchy_pdf (double x);
>   double test_chisq (void);
>   double test_chisq_pdf (double x);
>   double test_dirichlet (void);
>   double test_dirichlet_pdf (double x);
>   void test_dirichlet_moments (void);
> --- 35,64 ----
>   
>   void testMoments (double (*f) (void), const char *name,
>   		  double a, double b, double p);
> ! void testPDF (double (*f) (void), double (*pdf) (double),
> ! 	      double (*cdf) (double), const char *name);
>   void testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
> ! 		      double (*cdf) (unsigned int), const char *name);
>   
>   void test_shuffle (void);
>   void test_choose (void);
>   double test_beta (void);
>   double test_beta_pdf (double x);
> + double test_beta_cdf (double x);
>   double test_bernoulli (void);
>   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_cauchy (void);
>   double test_cauchy_pdf (double x);
> + double test_cauchy_cdf (double x);
>   double test_chisq (void);
>   double test_chisq_pdf (double x);
> + double test_chisq_cdf (double x);
>   double test_dirichlet (void);
>   double test_dirichlet_pdf (double x);
>   void test_dirichlet_moments (void);
> ***************
> *** 64,69 ****
> --- 70,76 ----
>   double test_erlang_pdf (double x);
>   double test_exponential (void);
>   double test_exponential_pdf (double x);
> + double test_exponential_cdf (double x);
>   double test_exppow0 (void);
>   double test_exppow0_pdf (double x);
>   double test_exppow1 (void);
> ***************
> *** 80,93 ****
> --- 87,103 ----
>   double test_flat_pdf (double x);
>   double test_gamma (void);
>   double test_gamma_pdf (double x);
> + double test_gamma_cdf (double x);
>   double test_gamma1 (void);
>   double test_gamma1_pdf (double x);
> + double test_gamma1_cdf (double x);
>   double test_gamma_int (void);
>   double test_gamma_int_pdf (double x);
>   double test_gamma_large (void);
>   double test_gamma_large_pdf (double x);
>   double test_gaussian (void);
>   double test_gaussian_pdf (double x);
> + double test_gaussian_cdf (double x);
>   double test_gaussian_ratio_method (void);
>   double test_gaussian_ratio_method_pdf (double x);
>   double test_gaussian_tail (void);
> ***************
> *** 116,123 ****
> --- 126,135 ----
>   double test_gumbel2_pdf (double x);
>   double test_geometric (void);
>   double test_geometric_pdf (unsigned int x);
> + double test_geometric_cdf (unsigned int x);
>   double test_geometric1 (void);
>   double test_geometric1_pdf (unsigned int x);
> + double test_geometric1_cdf (unsigned int x);
>   double test_hypergeometric1 (void);
>   double test_hypergeometric1_pdf (unsigned int x);
>   double test_hypergeometric2 (void);
> ***************
> *** 154,159 ****
> --- 166,172 ----
>   double test_levy_skew2b_pdf (double x);
>   double test_logistic (void);
>   double test_logistic_pdf (double x);
> + double test_logistic_cdf (double x);
>   double test_lognormal (void);
>   double test_lognormal_pdf (double x);
>   double test_logarithmic (void);
> ***************
> *** 169,174 ****
> --- 182,188 ----
>   double test_pascal_pdf (unsigned int n);
>   double test_pareto (void);
>   double test_pareto_pdf (double x);
> + double test_pareto_cdf (double x);
>   double test_poisson (void);
>   double test_poisson_pdf (unsigned int x);
>   double test_poisson_large (void);
> ***************
> *** 189,196 ****
> --- 203,212 ----
>   double test_rayleigh_tail_pdf (double x);
>   double test_tdist1 (void);
>   double test_tdist1_pdf (double x);
> + double test_tdist1_cdf (double x);
>   double test_tdist2 (void);
>   double test_tdist2_pdf (double x);
> + double test_tdist2_cdf (double x);
>   double test_laplace (void);
>   double test_laplace_pdf (double x);
>   double test_weibull (void);
> ***************
> *** 209,215 ****
>     r_global = gsl_rng_alloc (gsl_rng_default);
>   
>   #define FUNC(x)  test_ ## x,                     "test gsl_ran_" #x
> ! #define FUNC2(x) test_ ## x, test_ ## x ## _pdf, "test gsl_ran_" #x
>   
>     test_shuffle ();
>     test_choose ();
> --- 225,233 ----
>     r_global = gsl_rng_alloc (gsl_rng_default);
>   
>   #define FUNC(x)  test_ ## x,                     "test gsl_ran_" #x
> ! #define FUNC2(x) test_ ## x, test_ ## x ## _pdf, NULL, "test gsl_ran_" #x
> ! #define FUNC3(x) test_ ## x, test_ ## x ## _pdf, test_ ## x ## _cdf, \
> ! 		 "test gsl_ran_" #x
>   
>     test_shuffle ();
>     test_choose ();
> ***************
> *** 228,239 ****
>     test_dirichlet_moments ();
>     test_multinomial_moments ();
>   
> !   testPDF (FUNC2 (beta));
> !   testPDF (FUNC2 (cauchy));
> !   testPDF (FUNC2 (chisq));
>     testPDF (FUNC2 (dirichlet));
>     testPDF (FUNC2 (erlang));
> !   testPDF (FUNC2 (exponential));
>   
>     testPDF (FUNC2 (exppow0));
>     testPDF (FUNC2 (exppow1));
> --- 246,257 ----
>     test_dirichlet_moments ();
>     test_multinomial_moments ();
>   
> !   testPDF (FUNC3 (beta));	/* CDF available */
> !   testPDF (FUNC3 (cauchy));	/* CDF available */
> !   testPDF (FUNC3 (chisq));	/* CDF available */
>     testPDF (FUNC2 (dirichlet));
>     testPDF (FUNC2 (erlang));
> !   testPDF (FUNC3 (exponential));/* CDF available */
>   
>     testPDF (FUNC2 (exppow0));
>     testPDF (FUNC2 (exppow1));
> ***************
> *** 243,253 ****
>   
>     testPDF (FUNC2 (fdist));
>     testPDF (FUNC2 (flat));
> !   testPDF (FUNC2 (gamma));
> !   testPDF (FUNC2 (gamma1));
>     testPDF (FUNC2 (gamma_int));
>     testPDF (FUNC2 (gamma_large));
> !   testPDF (FUNC2 (gaussian));
>     testPDF (FUNC2 (gaussian_ratio_method));
>     testPDF (FUNC2 (ugaussian));
>     testPDF (FUNC2 (ugaussian_ratio_method));
> --- 261,271 ----
>   
>     testPDF (FUNC2 (fdist));
>     testPDF (FUNC2 (flat));
> !   testPDF (FUNC3 (gamma));	/* CDF available */
> !   testPDF (FUNC3 (gamma1));	/* CDF available */
>     testPDF (FUNC2 (gamma_int));
>     testPDF (FUNC2 (gamma_large));
> !   testPDF (FUNC3 (gaussian));  /* CDF available */
>     testPDF (FUNC2 (gaussian_ratio_method));
>     testPDF (FUNC2 (ugaussian));
>     testPDF (FUNC2 (ugaussian_ratio_method));
> ***************
> *** 274,286 ****
>     testPDF (FUNC2 (levy_skew2a));
>     testPDF (FUNC2 (levy_skew1b));
>     testPDF (FUNC2 (levy_skew2b));
> !   testPDF (FUNC2 (logistic));
>     testPDF (FUNC2 (lognormal));
> !   testPDF (FUNC2 (pareto));
>     testPDF (FUNC2 (rayleigh));
>     testPDF (FUNC2 (rayleigh_tail));
> !   testPDF (FUNC2 (tdist1));
> !   testPDF (FUNC2 (tdist2));
>     testPDF (FUNC2 (laplace));
>     testPDF (FUNC2 (weibull));
>     testPDF (FUNC2 (weibull1));
> --- 292,304 ----
>     testPDF (FUNC2 (levy_skew2a));
>     testPDF (FUNC2 (levy_skew1b));
>     testPDF (FUNC2 (levy_skew2b));
> !   testPDF (FUNC3 (logistic));	/* CDF available */
>     testPDF (FUNC2 (lognormal));
> !   testPDF (FUNC3 (pareto));	/* CDF available */
>     testPDF (FUNC2 (rayleigh));
>     testPDF (FUNC2 (rayleigh_tail));
> !   testPDF (FUNC3 (tdist1));	/* CDF available */
> !   testPDF (FUNC3 (tdist2));	/* CDF available */
>     testPDF (FUNC2 (laplace));
>     testPDF (FUNC2 (weibull));
>     testPDF (FUNC2 (weibull1));
> ***************
> *** 296,305 ****
>     testDiscretePDF (FUNC2 (poisson));
>     testDiscretePDF (FUNC2 (poisson_large));
>     testDiscretePDF (FUNC2 (bernoulli));
> !   testDiscretePDF (FUNC2 (binomial));
> !   testDiscretePDF (FUNC2 (binomial_large));
> !   testDiscretePDF (FUNC2 (geometric));
> !   testDiscretePDF (FUNC2 (geometric1));
>     testDiscretePDF (FUNC2 (hypergeometric1));
>     testDiscretePDF (FUNC2 (hypergeometric2));
>     testDiscretePDF (FUNC2 (hypergeometric3));
> --- 314,323 ----
>     testDiscretePDF (FUNC2 (poisson));
>     testDiscretePDF (FUNC2 (poisson_large));
>     testDiscretePDF (FUNC2 (bernoulli));
> !   testDiscretePDF (FUNC3 (binomial));  /* CDF available */
> !   testDiscretePDF (FUNC3 (binomial_large));  /* CDF available */
> !   testDiscretePDF (FUNC3 (geometric));  /* CDF available */
> !   testDiscretePDF (FUNC3 (geometric1));  /* CDF available */
>     testDiscretePDF (FUNC2 (hypergeometric1));
>     testDiscretePDF (FUNC2 (hypergeometric2));
>     testDiscretePDF (FUNC2 (hypergeometric3));
> ***************
> *** 434,444 ****
>   #define BINS 100
>   
>   void
> ! testPDF (double (*f) (void), double (*pdf) (double), const char *name)
>   {
>     double count[BINS], p[BINS];
> !   double a = -5.0, b = +5.0;
> !   double dx = (b - a) / BINS;
>     int i, j, status = 0, status_i = 0;
>   
>     for (i = 0; i < BINS; i++)
> --- 452,466 ----
>   #define BINS 100
>   
>   void
> ! testPDF (double (*f) (void),
> ! 	 double (*pdf) (double),
> ! 	 double (*cdf) (double),
> ! 	 const char *name)
>   {
>     double count[BINS], p[BINS];
> !   const double a = -5.0, b = +5.0;
> !   const double dx = (b - a) / BINS;
> !   double integral = 0;
>     int i, j, status = 0, status_i = 0;
>   
>     for (i = 0; i < BINS; i++)
> ***************
> *** 470,475 ****
> --- 492,515 ----
>   	sum += pdf (x + j * dx / STEPS);
>   
>         p[i] = 0.5 * (pdf (x) + 2 * sum + pdf (x + dx - 1e-7)) * dx / STEPS;
> +       integral += p[i];
> + 
> +       if (cdf != NULL)
> + 	/* leave this test in place as long as some _cdf() functions
> + 	   are missing */
> + 	{
> + 	  /* check whether the pdf integrates up to the cdf */
> + 	  double exact_integral = cdf (x + dx) - cdf (x);
> + 	  gsl_test_rel(exact_integral, p[i], 1e-3,
> + 		       "%s integral on [%g,%g)", name, x, x + dx);
> + 	}
> +     }
> + 
> +   if (cdf != NULL)
> +     {
> +       double exact_integral = cdf (b) - cdf (a);
> +       gsl_test_rel(exact_integral, integral, 1e-6,
> + 		   "%s integral on [%g,%g)", name, a, b);
>       }
>   
>     for (i = 0; i < BINS; i++)
> ***************
> *** 498,507 ****
>   
>   void
>   testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
> ! 		 const char *name)
>   {
>     double count[BINS], p[BINS];
>     unsigned int i;
>     int status = 0, status_i = 0;
>   
>     for (i = 0; i < BINS; i++)
> --- 538,548 ----
>   
>   void
>   testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
> ! 		 double (*cdf) (unsigned int), const char *name)
>   {
>     double count[BINS], p[BINS];
>     unsigned int i;
> +   double sum = 0;
>     int status = 0, status_i = 0;
>   
>     for (i = 0; i < BINS; i++)
> ***************
> *** 514,521 ****
>   	count[r]++;
>       }
>   
> !   for (i = 0; i < BINS; i++)
>       p[i] = pdf (i);
>   
>     for (i = 0; i < BINS; i++)
>       {
> --- 555,566 ----
>   	count[r]++;
>       }
>   
> !   for (i = 0; i < BINS; i++) {
>       p[i] = pdf (i);
> +     sum += p[i];
> +     if (cdf != NULL)
> +       gsl_test_rel(cdf(i), sum, 1e-6, "%s sum_0^%d", name, i);
> +   }
>   
>     for (i = 0; i < BINS; i++)
>       {
> ***************
> *** 553,558 ****
> --- 598,609 ----
>   {
>     return gsl_ran_beta_pdf (x, 2.0, 3.0);
>   }
> + double
> + test_beta_cdf (double x)
> + {
> +   return gsl_ran_beta_cdf (x, 2.0, 3.0);
> + }
> + 
>   
>   double
>   test_bernoulli (void)
> ***************
> *** 580,585 ****
> --- 631,642 ----
>   }
>   
>   double
> + test_binomial_cdf (unsigned int n)
> + {
> +   return gsl_ran_binomial_cdf (n, 0.3, 5);
> + }
> + 
> + double
>   test_binomial_large (void)
>   {
>     return gsl_ran_binomial (r_global, 0.3, 55);
> ***************
> *** 590,595 ****
> --- 647,658 ----
>   {
>     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_cauchy (void)
> ***************
> *** 604,609 ****
> --- 667,679 ----
>   }
>   
>   double
> + test_cauchy_cdf (double x)
> + {
> +   return gsl_ran_cauchy_cdf (x, 2.0);
> + }
> + 
> + 
> + double
>   test_chisq (void)
>   {
>     return gsl_ran_chisq (r_global, 13.0);
> ***************
> *** 616,621 ****
> --- 686,697 ----
>   }
>   
>   double
> + test_chisq_cdf (double x)
> + {
> +   return gsl_ran_chisq_cdf (x, 13.0);
> + }
> + 
> + double
>   test_dir2d (void)
>   {
>     double x = 0, y = 0, theta;
> ***************
> *** 911,916 ****
> --- 987,998 ----
>   }
>   
>   double
> + test_exponential_cdf (double x)
> + {
> +   return gsl_ran_exponential_cdf (x, 2.0);
> + }
> + 
> + double
>   test_exppow0 (void)
>   {
>     return gsl_ran_exppow (r_global, 3.7, 0.3);
> ***************
> *** 1008,1013 ****
> --- 1090,1101 ----
>   }
>   
>   double
> + test_gamma_cdf (double x)
> + {
> +   return gsl_ran_gamma_cdf (x, 2.5, 2.17);
> + }
> + 
> + double
>   test_gamma1 (void)
>   {
>     return gsl_ran_gamma (r_global, 1.0, 2.17);
> ***************
> *** 1019,1024 ****
> --- 1107,1117 ----
>     return gsl_ran_gamma_pdf (x, 1.0, 2.17);
>   }
>   
> + double
> + test_gamma1_cdf (double x)
> + {
> +   return gsl_ran_gamma_cdf (x, 1.0, 2.17);
> + }
>   
>   double
>   test_gamma_int (void)
> ***************
> *** 1057,1062 ****
> --- 1150,1161 ----
>   {
>     return gsl_ran_gaussian_pdf (x, 3.0);
>   }
> + double
> + 
> + test_gaussian_cdf (double x)
> + {
> +   return gsl_ran_gaussian_cdf (x, 3.0);
> + }
>   
>   double
>   test_gaussian_ratio_method (void)
> ***************
> *** 1232,1237 ****
> --- 1331,1342 ----
>   }
>   
>   double
> + test_geometric_cdf (unsigned int n)
> + {
> +   return gsl_ran_geometric_cdf (n, 0.5);
> + }
> + 
> + double
>   test_geometric1 (void)
>   {
>     return gsl_ran_geometric (r_global, 1.0);
> ***************
> *** 1244,1249 ****
> --- 1349,1360 ----
>   }
>   
>   double
> + test_geometric1_cdf (unsigned int n)
> + {
> +   return gsl_ran_geometric_cdf (n, 1.0);
> + }
> + 
> + double
>   test_hypergeometric1 (void)
>   {
>     return gsl_ran_hypergeometric (r_global, 5, 7, 4);
> ***************
> *** 1491,1496 ****
> --- 1602,1614 ----
>   }
>   
>   double
> + test_logistic_cdf (double x)
> + {
> +   return gsl_ran_logistic_cdf (x, 3.1);
> + }
> + 
> + 
> + double
>   test_logarithmic (void)
>   {
>     return gsl_ran_logarithmic (r_global, 0.4);
> ***************
> *** 1532,1538 ****
>   double
>   test_multinomial_pdf (unsigned int n_0)
>   {
> !   /* The margional distribution of just 1 variate  is binomial. */
>     size_t K = 2;
>     /* Test use of weights instead of probabilities */
>     double p[] = { 0.4, 1.6};
> --- 1650,1656 ----
>   double
>   test_multinomial_pdf (unsigned int n_0)
>   {
> !   /* The marginal distribution of just 1 variate is binomial. */
>     size_t K = 2;
>     /* Test use of weights instead of probabilities */
>     double p[] = { 0.4, 1.6};
> ***************
> *** 1603,1608 ****
> --- 1721,1733 ----
>   }
>   
>   double
> + test_pareto_cdf (double x)
> + {
> +   return gsl_ran_pareto_cdf (x, 1.9, 2.75);
> + }
> + 
> + 
> + double
>   test_rayleigh (void)
>   {
>     return gsl_ran_rayleigh (r_global, 1.9);
> ***************
> *** 1665,1670 ****
> --- 1790,1801 ----
>   }
>   
>   double
> + test_tdist1_cdf (double x)
> + {
> +   return gsl_ran_tdist_cdf (x, 1.75);
> + }
> + 
> + double
>   test_tdist2 (void)
>   {
>     return gsl_ran_tdist (r_global, 12.75);
> ***************
> *** 1674,1679 ****
> --- 1805,1816 ----
>   test_tdist2_pdf (double x)
>   {
>     return gsl_ran_tdist_pdf (x, 12.75);
> + }
> + 
> + double
> + test_tdist2_cdf (double x)
> + {
> +   return gsl_ran_tdist_cdf (x, 12.75);
>   }
>   
>   



More information about the Gsl-discuss mailing list