1 /* statistic_tests.cpp file
2 *
3 * Copyright Jens Maurer 2000, 2002
4 * Distributed under the Boost Software License, Version 1.0. (See
5 * accompanying file LICENSE_1_0.txt or copy at
6 * http://www.boost.org/LICENSE_1_0.txt)
7 *
8 * $Id$
9 *
10 * Revision history
11 */
12
13 #include <iostream>
14 #include <iomanip>
15 #include <string>
16 #include <functional>
17 #include <vector>
18 #include <set>
19 #include <algorithm>
20
21 #include <boost/cstdint.hpp>
22 #include <boost/random.hpp>
23
24 #include <boost/math/special_functions/gamma.hpp>
25
26 #include <boost/math/distributions/uniform.hpp>
27 #include <boost/math/distributions/chi_squared.hpp>
28 #include <boost/math/distributions/normal.hpp>
29 #include <boost/math/distributions/triangular.hpp>
30 #include <boost/math/distributions/cauchy.hpp>
31 #include <boost/math/distributions/gamma.hpp>
32 #include <boost/math/distributions/exponential.hpp>
33 #include <boost/math/distributions/lognormal.hpp>
34
35 #include "statistic_tests.hpp"
36 #include "integrate.hpp"
37
38 class test_environment;
39
40 class test_base
41 {
42 protected:
test_base(test_environment & env)43 explicit test_base(test_environment & env) : environment(env) { }
44 void check_(double val) const;
45 private:
46 test_environment & environment;
47 };
48
49 class equidistribution_test : test_base
50 {
51 public:
equidistribution_test(test_environment & env,unsigned int classes,unsigned int high_classes)52 equidistribution_test(test_environment & env, unsigned int classes,
53 unsigned int high_classes)
54 : test_base(env), classes(classes),
55 test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
56 { }
57
58 template<class RNG>
run(RNG & rng,int n1,int n2)59 void run(RNG & rng, int n1, int n2)
60 {
61 using namespace boost;
62 std::cout << "equidistribution: " << std::flush;
63 equidistribution_experiment equi(classes);
64 variate_generator<RNG&, uniform_smallint<> > uint_linear(rng, uniform_smallint<>(0, classes-1));
65 check_(run_experiment(test_distrib_chi_square,
66 experiment_generator(equi, uint_linear, n1), n2));
67 check_(run_experiment(test_distrib_chi_square,
68 experiment_generator(equi, uint_linear, n1), 2*n2));
69
70 std::cout << " 2D: " << std::flush;
71 equidistribution_2d_experiment equi_2d(classes);
72 unsigned int root = static_cast<unsigned int>(std::sqrt(double(classes)));
73 assert(root * root == classes);
74 variate_generator<RNG&, uniform_smallint<> > uint_square(rng, uniform_smallint<>(0, root-1));
75 check_(run_experiment(test_distrib_chi_square,
76 experiment_generator(equi_2d, uint_square, n1), n2));
77 check_(run_experiment(test_distrib_chi_square,
78 experiment_generator(equi_2d, uint_square, n1), 2*n2));
79 std::cout << std::endl;
80 }
81 private:
82 unsigned int classes;
83 distribution_experiment test_distrib_chi_square;
84 };
85
86 class ks_distribution_test : test_base
87 {
88 public:
ks_distribution_test(test_environment & env,unsigned int classes)89 ks_distribution_test(test_environment & env, unsigned int classes)
90 : test_base(env),
91 test_distrib_chi_square(kolmogorov_smirnov_probability(5000),
92 classes)
93 { }
94
95 template<class RNG>
run(RNG & rng,int n1,int n2)96 void run(RNG & rng, int n1, int n2)
97 {
98 boost::math::uniform ud(static_cast<double>((rng.min)()), static_cast<double>((rng.max)()));
99 run(rng, ud, n1, n2);
100 }
101 template<class RNG, class Dist>
run(RNG & rng,const Dist & dist,int n1,int n2)102 void run(RNG & rng, const Dist& dist, int n1, int n2)
103 {
104 using namespace boost;
105 std::cout << "KS: " << std::flush;
106 kolmogorov_experiment ks(n1);
107 check_(run_experiment(test_distrib_chi_square,
108 ks_experiment_generator(ks, rng, dist), n2));
109 check_(run_experiment(test_distrib_chi_square,
110 ks_experiment_generator(ks, rng, dist), 2*n2));
111 std::cout << std::endl;
112 }
113 private:
114 distribution_experiment test_distrib_chi_square;
115 };
116
117 class runs_test : test_base
118 {
119 public:
runs_test(test_environment & env,unsigned int classes,unsigned int high_classes)120 runs_test(test_environment & env, unsigned int classes,
121 unsigned int high_classes)
122 : test_base(env), classes(classes),
123 test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
124 { }
125
126 template<class RNG>
run(RNG & rng,int n1,int n2)127 void run(RNG & rng, int n1, int n2)
128 {
129 using namespace boost;
130 std::cout << "runs: up: " << std::flush;
131 runs_experiment<true> r_up(classes);
132
133 check_(run_experiment(test_distrib_chi_square,
134 experiment_generator(r_up, rng, n1), n2));
135 check_(run_experiment(test_distrib_chi_square,
136 experiment_generator(r_up, rng, n1), 2*n2));
137
138 std::cout << " down: " << std::flush;
139 runs_experiment<false> r_down(classes);
140 check_(run_experiment(test_distrib_chi_square,
141 experiment_generator(r_down, rng, n1), n2));
142 check_(run_experiment(test_distrib_chi_square,
143 experiment_generator(r_down, rng, n1), 2*n2));
144
145 std::cout << std::endl;
146 }
147 private:
148 unsigned int classes;
149 distribution_experiment test_distrib_chi_square;
150 };
151
152 class gap_test : test_base
153 {
154 public:
gap_test(test_environment & env,unsigned int classes,unsigned int high_classes)155 gap_test(test_environment & env, unsigned int classes,
156 unsigned int high_classes)
157 : test_base(env), classes(classes),
158 test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
159 { }
160
161 template<class RNG>
run(RNG & rng,int n1,int n2)162 void run(RNG & rng, int n1, int n2)
163 {
164 boost::math::uniform ud(
165 static_cast<double>((rng.min)()),
166 static_cast<double>((rng.max)()) +
167 (std::numeric_limits<typename RNG::result_type>::is_integer? 0.0 : 1.0));
168 run(rng, ud, n1, n2);
169 }
170
171 template<class RNG, class Dist>
run(RNG & rng,const Dist & dist,int n1,int n2)172 void run(RNG & rng, const Dist& dist, int n1, int n2)
173 {
174 using namespace boost;
175 std::cout << "gaps: " << std::flush;
176 gap_experiment gap(classes, dist, 0.2, 0.8);
177
178 check_(run_experiment(test_distrib_chi_square,
179 experiment_generator(gap, rng, n1), n2));
180 check_(run_experiment(test_distrib_chi_square,
181 experiment_generator(gap, rng, n1), 2*n2));
182
183 std::cout << std::endl;
184 }
185 private:
186 unsigned int classes;
187 distribution_experiment test_distrib_chi_square;
188 };
189
190 class poker_test : test_base
191 {
192 public:
poker_test(test_environment & env,unsigned int classes,unsigned int high_classes)193 poker_test(test_environment & env, unsigned int classes,
194 unsigned int high_classes)
195 : test_base(env), classes(classes),
196 test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
197 { }
198
199 template<class RNG>
run(RNG & rng,int n1,int n2)200 void run(RNG & rng, int n1, int n2)
201 {
202 using namespace boost;
203 std::cout << "poker: " << std::flush;
204 poker_experiment poker(8, classes);
205 variate_generator<RNG&, uniform_smallint<> > usmall(rng, uniform_smallint<>(0, 7));
206 check_(run_experiment(test_distrib_chi_square,
207 experiment_generator(poker, usmall, n1), n2));
208 check_(run_experiment(test_distrib_chi_square,
209 experiment_generator(poker, usmall, n1), 2*n2));
210 std::cout << std::endl;
211 }
212 private:
213 unsigned int classes;
214 distribution_experiment test_distrib_chi_square;
215 };
216
217 class coupon_collector_test : test_base
218 {
219 public:
coupon_collector_test(test_environment & env,unsigned int classes,unsigned int high_classes)220 coupon_collector_test(test_environment & env, unsigned int classes,
221 unsigned int high_classes)
222 : test_base(env), classes(classes),
223 test_distrib_chi_square(boost::math::chi_squared(classes-1), high_classes)
224 { }
225
226 template<class RNG>
run(RNG & rng,int n1,int n2)227 void run(RNG & rng, int n1, int n2)
228 {
229 using namespace boost;
230 std::cout << "coupon collector: " << std::flush;
231 coupon_collector_experiment coupon(5, classes);
232
233 variate_generator<RNG&, uniform_smallint<> > usmall(rng, uniform_smallint<>(0, 4));
234 check_(run_experiment(test_distrib_chi_square,
235 experiment_generator(coupon, usmall, n1), n2));
236 check_(run_experiment(test_distrib_chi_square,
237 experiment_generator(coupon, usmall, n1), 2*n2));
238 std::cout << std::endl;
239 }
240 private:
241 unsigned int classes;
242 distribution_experiment test_distrib_chi_square;
243 };
244
245 class permutation_test : test_base
246 {
247 public:
permutation_test(test_environment & env,unsigned int classes,unsigned int high_classes)248 permutation_test(test_environment & env, unsigned int classes,
249 unsigned int high_classes)
250 : test_base(env), classes(classes),
251 test_distrib_chi_square(boost::math::chi_squared(fac<int>(classes)-1),
252 high_classes)
253 { }
254
255 template<class RNG>
run(RNG & rng,int n1,int n2)256 void run(RNG & rng, int n1, int n2)
257 {
258 using namespace boost;
259 std::cout << "permutation: " << std::flush;
260 permutation_experiment perm(classes);
261
262 // generator_reference_t<RNG> gen_ref(rng);
263 RNG& gen_ref(rng);
264 check_(run_experiment(test_distrib_chi_square,
265 experiment_generator(perm, gen_ref, n1), n2));
266 check_(run_experiment(test_distrib_chi_square,
267 experiment_generator(perm, gen_ref, n1), 2*n2));
268 std::cout << std::endl;
269 }
270 private:
271 unsigned int classes;
272 distribution_experiment test_distrib_chi_square;
273 };
274
275 class maximum_test : test_base
276 {
277 public:
maximum_test(test_environment & env,unsigned int high_classes)278 maximum_test(test_environment & env, unsigned int high_classes)
279 : test_base(env),
280 test_distrib_chi_square(kolmogorov_smirnov_probability(1000),
281 high_classes)
282 { }
283
284 template<class RNG>
run(RNG & rng,int n1,int n2)285 void run(RNG & rng, int n1, int n2)
286 {
287 using namespace boost;
288 std::cout << "maximum-of-t: " << std::flush;
289 maximum_experiment<RNG> mx(rng, n1, 5);
290 check_(run_experiment(test_distrib_chi_square, mx, n2));
291 check_(run_experiment(test_distrib_chi_square, mx, 2*n2));
292 std::cout << std::endl;
293 }
294 private:
295 distribution_experiment test_distrib_chi_square;
296 };
297
298 class birthday_test : test_base
299 {
300 public:
birthday_test(test_environment & env,unsigned int high_classes)301 birthday_test(test_environment & env, unsigned int high_classes)
302 : test_base(env),
303 test_distrib_chi_square(boost::math::chi_squared(4-1), high_classes)
304 { }
305
306 template<class RNG>
run(RNG & rng,int n1,int n2)307 void run(RNG & rng, int n1, int n2)
308 {
309 using namespace boost;
310 std::cout << "birthday spacing: " << std::flush;
311 boost::variate_generator<RNG&, boost::uniform_int<> > uni(rng, boost::uniform_int<>(0, (1<<25)-1));
312 birthday_spacing_experiment bsp(4, 512, (1<<25));
313 check_(run_experiment(test_distrib_chi_square,
314 experiment_generator(bsp, uni, n1), n2));
315 check_(run_experiment(test_distrib_chi_square,
316 experiment_generator(bsp, uni, n1), 2*n2));
317 std::cout << std::endl;
318 }
319 private:
320 distribution_experiment test_distrib_chi_square;
321 };
322
323 #ifdef BOOST_MSVC
324 #pragma warning(disable:4355)
325 #endif
326
327 class test_environment
328 {
329 public:
330 static const int classes = 20;
test_environment(double confid)331 explicit test_environment(double confid)
332 : confidence(confid),
333 confidence_chi_square_quantil(quantile(boost::math::chi_squared(classes-1), confidence)),
334 test_distrib_chi_square6(boost::math::chi_squared(7-1), classes),
335 ksdist_test(*this, classes),
336 equi_test(*this, 100, classes),
337 rns_test(*this, 7, classes),
338 gp_test(*this, 7, classes),
339 pk_test(*this, 5, classes),
340 cpn_test(*this, 15, classes),
341 perm_test(*this, 5, classes),
342 max_test(*this, classes),
343 bday_test(*this, classes)
344 {
345 std::cout << "Confidence level: " << confid
346 << "; 1-alpha = " << (1-confid)
347 << "; chi_square(" << (classes-1)
348 << ", " << confidence_chi_square_quantil
349 << ") = "
350 << cdf(boost::math::chi_squared(classes-1), confidence_chi_square_quantil)
351 << std::endl;
352 }
353
check_confidence(double val,double chi_square_conf) const354 bool check_confidence(double val, double chi_square_conf) const
355 {
356 std::cout << val;
357 bool result = (val <= chi_square_conf);
358 if(!result) {
359 std::cout << "* [";
360 double prob = (val > 10*chi_square_conf ? 1 :
361 cdf(boost::math::chi_squared(classes-1), val));
362 std::cout << (1-prob) << "]";
363 }
364 std::cout << " " << std::flush;
365 return result;
366 }
367
check_(double chi_square_value) const368 bool check_(double chi_square_value) const
369 {
370 return check_confidence(chi_square_value, confidence_chi_square_quantil);
371 }
372
373 template<class RNG>
run_test(const std::string & name)374 void run_test(const std::string & name)
375 {
376 using namespace boost;
377
378 std::cout << "Running tests on " << name << std::endl;
379
380 RNG rng(1234567);
381
382 ksdist_test.run(rng, 5000, 250);
383 equi_test.run(rng, 5000, 250);
384 rns_test.run(rng, 100000, 250);
385 gp_test.run(rng, 10000, 250);
386 pk_test.run(rng, 5000, 250);
387 cpn_test.run(rng, 500, 250);
388 perm_test.run(rng, 1200, 250);
389 max_test.run(rng, 1000, 250);
390 bday_test.run(rng, 1000, 150);
391
392 std::cout << std::endl;
393 }
394
395 template<class RNG, class Dist, class ExpectedDist>
run_test(const std::string & name,const Dist & dist,const ExpectedDist & expected_dist)396 void run_test(const std::string & name, const Dist & dist, const ExpectedDist & expected_dist)
397 {
398 using namespace boost;
399
400 std::cout << "Running tests on " << name << std::endl;
401
402 RNG rng;
403 variate_generator<RNG&, Dist> vgen(rng, dist);
404
405 ksdist_test.run(vgen, expected_dist, 5000, 250);
406 rns_test.run(vgen, 100000, 250);
407 gp_test.run(vgen, expected_dist, 10000, 250);
408 perm_test.run(vgen, 1200, 250);
409
410 std::cout << std::endl;
411 }
412
413 private:
414 double confidence;
415 double confidence_chi_square_quantil;
416 distribution_experiment test_distrib_chi_square6;
417 ks_distribution_test ksdist_test;
418 equidistribution_test equi_test;
419 runs_test rns_test;
420 gap_test gp_test;
421 poker_test pk_test;
422 coupon_collector_test cpn_test;
423 permutation_test perm_test;
424 maximum_test max_test;
425 birthday_test bday_test;
426 };
427
check_(double val) const428 void test_base::check_(double val) const
429 {
430 environment.check_(val);
431 }
432
433 class program_args
434 {
435 public:
program_args(int argc,char ** argv)436 program_args(int argc, char** argv)
437 {
438 if(argc > 0) {
439 names.insert(argv + 1, argv + argc);
440 }
441 }
check_(const std::string & test_name) const442 bool check_(const std::string & test_name) const
443 {
444 return(names.empty() || names.find(test_name) != names.end());
445 }
446 private:
447 std::set<std::string> names;
448 };
449
main(int argc,char * argv[])450 int main(int argc, char* argv[])
451 {
452 program_args args(argc, argv);
453 test_environment env(0.99);
454
455 #define TEST(name) \
456 if(args.check_(#name)) \
457 env.run_test<boost::name>(#name)
458
459 TEST(minstd_rand0);
460 TEST(minstd_rand);
461 TEST(rand48);
462 TEST(ecuyer1988);
463 TEST(kreutzer1986);
464 TEST(taus88);
465 TEST(hellekalek1995);
466 TEST(mt11213b);
467 TEST(mt19937);
468 TEST(lagged_fibonacci607);
469 TEST(lagged_fibonacci1279);
470 TEST(lagged_fibonacci2281);
471 TEST(lagged_fibonacci3217);
472 TEST(lagged_fibonacci4423);
473 TEST(lagged_fibonacci9689);
474 TEST(lagged_fibonacci19937);
475 TEST(lagged_fibonacci23209);
476 TEST(lagged_fibonacci44497);
477 TEST(ranlux3);
478 TEST(ranlux4);
479
480 #if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
481 TEST(ranlux64_3);
482 TEST(ranlux64_4);
483 #endif
484
485 TEST(ranlux3_01);
486 TEST(ranlux4_01);
487 TEST(ranlux64_3_01);
488 TEST(ranlux64_4_01);
489
490 if(args.check_("normal"))
491 env.run_test<boost::mt19937>("normal", boost::normal_distribution<>(), boost::math::normal());
492 if(args.check_("triangle"))
493 env.run_test<boost::mt19937>("triangle", boost::triangle_distribution<>(0, 1, 3), boost::math::triangular(0, 1, 3));
494 if(args.check_("cauchy"))
495 env.run_test<boost::mt19937>("cauchy", boost::cauchy_distribution<>(), boost::math::cauchy());
496 if(args.check_("gamma"))
497 env.run_test<boost::mt19937>("gamma", boost::gamma_distribution<>(1), boost::math::gamma_distribution<>(1));
498 if(args.check_("exponential"))
499 env.run_test<boost::mt19937>("exponential", boost::exponential_distribution<>(), boost::math::exponential());
500 if(args.check_("lognormal"))
501 env.run_test<boost::mt19937>("lognormal", boost::lognormal_distribution<>(1, 1),
502 boost::math::lognormal(std::log(1.0/std::sqrt(2.0)), std::sqrt(std::log(2.0))));
503 }
504