Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r60335 - trunk/libs/random/test
From: steven_at_[hidden]
Date: 2010-03-07 22:46:48


Author: steven_watanabe
Date: 2010-03-07 22:46:47 EST (Sun, 07 Mar 2010)
New Revision: 60335
URL: http://svn.boost.org/trac/boost/changeset/60335

Log:
Add distributions to statistic_tests
Text files modified:
   trunk/libs/random/test/integrate.hpp | 8
   trunk/libs/random/test/statistic_tests.cpp | 162 ++++++++++++++++-----------------------
   trunk/libs/random/test/statistic_tests.hpp | 101 ++++++++++++++++++++----
   3 files changed, 153 insertions(+), 118 deletions(-)

Modified: trunk/libs/random/test/integrate.hpp
==============================================================================
--- trunk/libs/random/test/integrate.hpp (original)
+++ trunk/libs/random/test/integrate.hpp 2010-03-07 22:46:47 EST (Sun, 07 Mar 2010)
@@ -43,11 +43,11 @@
 }
 
 // compute b so that f(b) = y; assume f is monotone increasing
-template<class UnaryFunction>
-inline typename UnaryFunction::argument_type
+template<class UnaryFunction, class T>
+inline T
 invert_monotone_inc(UnaryFunction f, typename UnaryFunction::result_type y,
- typename UnaryFunction::argument_type lower = -1,
- typename UnaryFunction::argument_type upper = 1)
+ T lower = -1,
+ T upper = 1)
 {
   while(upper-lower > 1e-6) {
     double middle = (upper+lower)/2;

Modified: trunk/libs/random/test/statistic_tests.cpp
==============================================================================
--- trunk/libs/random/test/statistic_tests.cpp (original)
+++ trunk/libs/random/test/statistic_tests.cpp 2010-03-07 22:46:47 EST (Sun, 07 Mar 2010)
@@ -23,64 +23,18 @@
 
 #include <boost/math/special_functions/gamma.hpp>
 
+#include <boost/math/distributions/uniform.hpp>
+#include <boost/math/distributions/chi_squared.hpp>
+#include <boost/math/distributions/normal.hpp>
+#include <boost/math/distributions/triangular.hpp>
+#include <boost/math/distributions/cauchy.hpp>
+#include <boost/math/distributions/gamma.hpp>
+#include <boost/math/distributions/exponential.hpp>
+#include <boost/math/distributions/lognormal.hpp>
+
 #include "statistic_tests.hpp"
 #include "integrate.hpp"
 
-
-double normal_density(double x)
-{
- const double pi = 3.14159265358979323846;
- return 1/std::sqrt(2*pi) * std::exp(-x*x/2);
-}
-
-class chi_square_density : public std::unary_function<double, double>
-{
-public:
- chi_square_density(int freedom)
- : _exponent( static_cast<double>(freedom)/2-1 ),
- _factor(1/(std::pow(2, _exponent+1) * std::exp(boost::math::lgamma(_exponent+1))))
- { }
-
- double operator()(double x)
- {
- return _factor*std::pow(x, _exponent)*std::exp(-x/2);
- }
-private:
- double _exponent, _factor;
-};
-
-// computes F(x) or F(y) - F(x)
-class chi_square_probability : public distribution_function<double>
-{
-public:
- chi_square_probability(int freedom) : dens(freedom) {}
- double operator()(double x) { return operator()(0, x); }
- double operator()(double x, double y)
- { return trapezoid(dens, x, y, 1000); }
-private:
- chi_square_density dens;
-};
-
-class uniform_distribution : public distribution_function<double>
-{
-public:
- uniform_distribution(double from, double to) : from(from), to(to)
- { assert(from < to); }
- double operator()(double x)
- {
- if(x < from)
- return 0;
- else if(x > to)
- return 1;
- else
- return (x-from)/(to-from);
- }
- double operator()(double x, double delta)
- { return operator()(x+delta) - operator()(x); }
-private:
- double from, to;
-};
-
 class test_environment;
 
 class test_base
@@ -98,7 +52,7 @@
   equidistribution_test(test_environment & env, unsigned int classes,
                         unsigned int high_classes)
     : test_base(env), classes(classes),
- test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
+ test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
   { }
 
   template<class RNG>
@@ -129,10 +83,10 @@
   distribution_experiment test_distrib_chi_square;
 };
 
-class ks_equidistribution_test : test_base
+class ks_distribution_test : test_base
 {
 public:
- ks_equidistribution_test(test_environment & env, unsigned int classes)
+ ks_distribution_test(test_environment & env, unsigned int classes)
     : test_base(env),
       test_distrib_chi_square(kolmogorov_smirnov_probability(5000),
                               classes)
@@ -141,16 +95,19 @@
   template<class RNG>
   void run(RNG & rng, int n1, int n2)
   {
+ boost::math::uniform ud(static_cast<double>((rng.min)()), static_cast<double>((rng.max)()));
+ run(rng, ud, n1, n2);
+ }
+ template<class RNG, class Dist>
+ void run(RNG & rng, const Dist& dist, int n1, int n2)
+ {
     using namespace boost;
     std::cout << "KS: " << std::flush;
- // generator_reference_t<RNG> gen_ref(rng);
- RNG& gen_ref(rng);
     kolmogorov_experiment ks(n1);
- uniform_distribution ud(static_cast<double>((rng.min)()), static_cast<double>((rng.max)()));
     check(run_experiment(test_distrib_chi_square,
- ks_experiment_generator(ks, gen_ref, ud), n2));
+ ks_experiment_generator(ks, rng, dist), n2));
     check(run_experiment(test_distrib_chi_square,
- ks_experiment_generator(ks, gen_ref, ud), 2*n2));
+ ks_experiment_generator(ks, rng, dist), 2*n2));
     std::cout << std::endl;
   }
 private:
@@ -163,7 +120,7 @@
   runs_test(test_environment & env, unsigned int classes,
             unsigned int high_classes)
     : test_base(env), classes(classes),
- test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
+ test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
   { }
 
   template<class RNG>
@@ -199,7 +156,7 @@
   gap_test(test_environment & env, unsigned int classes,
             unsigned int high_classes)
     : test_base(env), classes(classes),
- test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
+ test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
   { }
 
   template<class RNG>
@@ -228,7 +185,7 @@
   poker_test(test_environment & env, unsigned int classes,
              unsigned int high_classes)
     : test_base(env), classes(classes),
- test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
+ test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
   { }
 
   template<class RNG>
@@ -255,7 +212,7 @@
   coupon_collector_test(test_environment & env, unsigned int classes,
                         unsigned int high_classes)
     : test_base(env), classes(classes),
- test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
+ test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
   { }
 
   template<class RNG>
@@ -283,7 +240,7 @@
   permutation_test(test_environment & env, unsigned int classes,
                    unsigned int high_classes)
     : test_base(env), classes(classes),
- test_distrib_chi_square(chi_square_probability(fac<int>(classes)-1),
+ test_distrib_chi_square(boost::math::chi_squared(fac<int>(classes)-1),
                               high_classes)
   { }
 
@@ -335,7 +292,7 @@
 public:
   birthday_test(test_environment & env, unsigned int high_classes)
     : test_base(env),
- test_distrib_chi_square(chi_square_probability(4-1), high_classes)
+ test_distrib_chi_square(boost::math::chi_squared(4-1), high_classes)
   { }
 
   template<class RNG>
@@ -365,9 +322,9 @@
   static const int classes = 20;
   explicit test_environment(double confid)
     : confidence(confid),
- confidence_chi_square_quantil(quantil(chi_square_density(classes-1), 0, confidence, 1e-4)),
- test_distrib_chi_square6(chi_square_probability(7-1), classes),
- ksequi_test(*this, classes),
+ confidence_chi_square_quantil(quantile(boost::math::chi_squared(classes-1), confidence)),
+ test_distrib_chi_square6(boost::math::chi_squared(7-1), classes),
+ ksdist_test(*this, classes),
       equi_test(*this, 100, classes),
       rns_test(*this, 7, classes),
       gp_test(*this, 7, classes),
@@ -382,7 +339,7 @@
               << "; chi_square(" << (classes-1)
               << ", " << confidence_chi_square_quantil
               << ") = "
- << chi_square_probability(classes-1)(0, confidence_chi_square_quantil)
+ << cdf(boost::math::chi_squared(classes-1), confidence_chi_square_quantil)
               << std::endl;
   }
   
@@ -393,7 +350,7 @@
     if(!result) {
       std::cout << "* [";
       double prob = (val > 10*chi_square_conf ? 1 :
- chi_square_probability(classes-1)(0, val));
+ cdf(boost::math::chi_squared(classes-1), val));
       std::cout << (1-prob) << "]";
     }
     std::cout << " " << std::flush;
@@ -413,10 +370,8 @@
     std::cout << "Running tests on " << name << std::endl;
 
     RNG rng(1234567);
- typedef boost::uniform_01<RNG> UGen;
 
-#if 1
- ksequi_test.run(rng, 5000, 250);
+ ksdist_test.run(rng, 5000, 250);
     equi_test.run(rng, 5000, 250);
     rns_test.run(rng, 100000, 250);
     gp_test.run(rng, 10000, 250);
@@ -424,17 +379,33 @@
     cpn_test.run(rng, 500, 250);
     perm_test.run(rng, 1200, 250);
     max_test.run(rng, 1000, 250);
-#endif
     bday_test.run(rng, 1000, 150);
 
     std::cout << std::endl;
   }
 
+ template<class RNG, class Dist, class ExpectedDist>
+ void run_test(const std::string & name, const Dist & dist, const ExpectedDist & expected_dist)
+ {
+ using namespace boost;
+
+ std::cout << "Running tests on " << name << std::endl;
+
+ RNG rng;
+ variate_generator<RNG&, Dist> vgen(rng, dist);
+
+ ksdist_test.run(vgen, expected_dist, 5000, 250);
+ rns_test.run(vgen, 100000, 250);
+ perm_test.run(vgen, 1200, 250);
+
+ std::cout << std::endl;
+ }
+
 private:
   double confidence;
   double confidence_chi_square_quantil;
   distribution_experiment test_distrib_chi_square6;
- ks_equidistribution_test ksequi_test;
+ ks_distribution_test ksdist_test;
   equidistribution_test equi_test;
   runs_test rns_test;
   gap_test gp_test;
@@ -450,23 +421,6 @@
   environment.check(val);
 }
 
-void print_ks_table()
-{
- std::cout.setf(std::ios::fixed);
- std::cout.precision(5);
- static const double all_p[] = { 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99 };
- for(int n = 0; n <= 10000; (n < 55 ? ++n : n *= 10)) {
- std::cout << std::setw(4) << n << " ";
- for(unsigned int i = 0; i < sizeof(all_p)/sizeof(all_p[0]); ++i) {
- std::cout << std::setw(8)
- << (n == 0 ? all_p[i] :
- invert_monotone_inc(kolmogorov_smirnov_probability(n), all_p[i], 0, 10))
- << " ";
- }
- std::cout << std::endl;
- }
-}
-
 class program_args
 {
 public:
@@ -476,7 +430,7 @@
       names.insert(argv + 1, argv + argc);
     }
   }
- bool check(std::string test_name) const
+ bool check(const std::string & test_name) const
   {
     return(names.empty() || names.find(test_name) != names.end());
   }
@@ -523,4 +477,18 @@
   TEST(ranlux4_01);
   TEST(ranlux64_3_01);
   TEST(ranlux64_4_01);
+
+ if(args.check("normal"))
+ env.run_test<boost::mt19937>("normal", boost::normal_distribution<>(), boost::math::normal());
+ if(args.check("triangle"))
+ env.run_test<boost::mt19937>("triangle", boost::triangle_distribution<>(0, 1, 3), boost::math::triangular(0, 1, 3));
+ if(args.check("cauchy"))
+ env.run_test<boost::mt19937>("cauchy", boost::cauchy_distribution<>(), boost::math::cauchy());
+ if(args.check("gamma"))
+ env.run_test<boost::mt19937>("gamma", boost::gamma_distribution<>(1), boost::math::gamma_distribution<>(1));
+ if(args.check("exponential"))
+ env.run_test<boost::mt19937>("exponential", boost::exponential_distribution<>(), boost::math::exponential());
+ if(args.check("lognormal"))
+ env.run_test<boost::mt19937>("lognormal", boost::lognormal_distribution<>(1, 1),
+ boost::math::lognormal(std::log(1.0/std::sqrt(2.0)), std::sqrt(std::log(2.0))));
 }

Modified: trunk/libs/random/test/statistic_tests.hpp
==============================================================================
--- trunk/libs/random/test/statistic_tests.hpp (original)
+++ trunk/libs/random/test/statistic_tests.hpp 2010-03-07 22:46:47 EST (Sun, 07 Mar 2010)
@@ -21,7 +21,9 @@
 
 #include <boost/random.hpp>
 #include <boost/config.hpp>
+#include <boost/bind.hpp>
 
+#include "integrate.hpp"
 
 #if defined(BOOST_MSVC) && BOOST_MSVC <= 1300
 namespace std
@@ -119,12 +121,12 @@
 class distribution_experiment : public equidistribution_experiment
 {
 public:
- template<class UnaryFunction>
- distribution_experiment(UnaryFunction probability , unsigned int classes)
+ template<class Distribution>
+ distribution_experiment(Distribution dist , unsigned int classes)
     : equidistribution_experiment(classes), limit(classes)
   {
     for(unsigned int i = 0; i < classes-1; ++i)
- limit[i] = invert_monotone_inc(probability, (i+1)*0.05, 0, 1000);
+ limit[i] = quantile(dist, (i+1)*0.05);
     limit[classes-1] = std::numeric_limits<double>::infinity();
     if(limit[classes-1] < (std::numeric_limits<double>::max)())
       limit[classes-1] = (std::numeric_limits<double>::max)();
@@ -150,6 +152,36 @@
   limits_type limit;
 };
 
+template<bool up, bool is_float>
+struct runs_direction_helper
+{
+ template<class T>
+ static T init(T)
+ {
+ return (std::numeric_limits<T>::max)();
+ }
+};
+
+template<>
+struct runs_direction_helper<true, true>
+{
+ template<class T>
+ static T init(T)
+ {
+ return -(std::numeric_limits<T>::max)();
+ }
+};
+
+template<>
+struct runs_direction_helper<true, false>
+{
+ template<class T>
+ static T init(T)
+ {
+ return (std::numeric_limits<T>::min)();
+ }
+};
+
 // runs-up/runs-down experiment
 template<bool up>
 class runs_experiment : public experiment_base
@@ -157,11 +189,15 @@
 public:
   explicit runs_experiment(unsigned int classes) : experiment_base(classes) { }
   
- template<class UniformRandomNumberGenerator, class Counter>
- void run(UniformRandomNumberGenerator & f, Counter & count, int n) const
+ template<class NumberGenerator, class Counter>
+ void run(NumberGenerator & f, Counter & count, int n) const
   {
- typedef typename UniformRandomNumberGenerator::result_type result_type;
- result_type init = (up ? (f.min)() : (f.max)());
+ typedef typename NumberGenerator::result_type result_type;
+ result_type init =
+ runs_direction_helper<
+ up,
+ !std::numeric_limits<result_type>::is_integer
+ >::init(result_type());
     result_type previous = init;
     unsigned int length = 0;
     for(int i = 0; i < n; ++i) {
@@ -423,7 +459,7 @@
       n_n = std::pow(static_cast<double>(n), n);
   }
   
- double operator()(double t) const
+ double cdf(double t) const
   {
     if(approx) {
       return 1-std::exp(-2*t*t)*(1-2.0/3.0*t/sqrt_n);
@@ -436,8 +472,8 @@
       return 1 - t/n_n * sum;
     }
   }
- double operator()(double t1, double t2) const
- { return operator()(t2) - operator()(t1); }
+ //double operator()(double t1, double t2) const
+ //{ return operator()(t2) - operator()(t1); }
 
 private:
   bool approx;
@@ -446,6 +482,16 @@
   double n_n;
 };
 
+inline double cdf(const kolmogorov_smirnov_probability& dist, double val)
+{
+ return dist.cdf(val);
+}
+
+inline double quantile(const kolmogorov_smirnov_probability& dist, double val)
+{
+ return invert_monotone_inc(boost::bind(&cdf, dist, _1), val, 0.0, 1000.0);
+}
+
 /*
  * Experiments for generators with continuous distribution functions
  */
@@ -462,7 +508,7 @@
     std::vector<int> c(m,0);
     for(int i = 0; i < n; ++i) {
       double val = static_cast<double>(gen());
- double y = distrib(val);
+ double y = cdf(distrib, val);
       int k = static_cast<int>(std::floor(m*y));
       if(k >= m)
         --k; // should not happen
@@ -486,13 +532,24 @@
   }
   double probability(double x) const
   {
- return ksp(x);
+ return cdf(ksp, x);
   }
 private:
   int n;
   kolmogorov_smirnov_probability ksp;
 };
 
+struct power_distribution
+{
+ power_distribution(double t) : t(t) {}
+ double t;
+};
+
+double cdf(const power_distribution& dist, double val)
+{
+ return std::pow(val, dist.t);
+}
+
 // maximum-of-t test (KS-based)
 template<class UniformRandomNumberGenerator>
 class maximum_experiment
@@ -504,9 +561,8 @@
 
   double operator()() const
   {
- double res = ke.run(generator(f, t),
- std::bind2nd(std::ptr_fun(static_cast<double (*)(double, double)>(&std::pow)), t));
- return res;
+ generator gen(f, t);
+ return ke.run(gen, power_distribution(t));
   }
 
 private:
@@ -590,6 +646,17 @@
                                        experiment));
 }
 
+// chi_square test
+template<class Experiment, class Generator>
+double run_experiment(const Experiment & experiment, const Generator & gen, int n)
+{
+ generic_counter<std::vector<int> > v(experiment.classes());
+ experiment.run(gen, v, n);
+ return chi_square_value(v.begin(), v.end(),
+ std::bind1st(std::mem_fun_ref(&Experiment::probability),
+ experiment));
+}
+
 // number generator with experiment results (for nesting)
 template<class Experiment, class Generator>
 class experiment_generator_t
@@ -597,7 +664,7 @@
 public:
   experiment_generator_t(const Experiment & exper, Generator & gen, int n)
     : experiment(exper), generator(gen), n(n) { }
- double operator()() { return run_experiment(experiment, generator, n); }
+ double operator()() const { return run_experiment(experiment, generator, n); }
 private:
   const Experiment & experiment;
   Generator & generator;
@@ -619,7 +686,7 @@
   ks_experiment_generator_t(const Experiment & exper, Generator & gen,
                             const Distribution & distrib)
     : experiment(exper), generator(gen), distribution(distrib) { }
- double operator()() { return experiment.run(generator, distribution); }
+ double operator()() const { return experiment.run(generator, distribution); }
 private:
   const Experiment & experiment;
   Generator & generator;


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk