[29] | 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: statistic_tests.cpp,v 1.11 2004/07/27 03:43:34 dgregor Exp $ |
---|
| 9 | * |
---|
| 10 | * Revision history |
---|
| 11 | */ |
---|
| 12 | |
---|
| 13 | /* |
---|
| 14 | * NOTE: This is not part of the official boost submission. It exists |
---|
| 15 | * only as a collection of ideas. |
---|
| 16 | */ |
---|
| 17 | |
---|
| 18 | #include <iostream> |
---|
| 19 | #include <iomanip> |
---|
| 20 | #include <string> |
---|
| 21 | #include <functional> |
---|
| 22 | #include <math.h> // lgamma is not in namespace std |
---|
| 23 | #include <vector> |
---|
| 24 | #include <algorithm> |
---|
| 25 | |
---|
| 26 | #include <boost/cstdint.hpp> |
---|
| 27 | #include <boost/random.hpp> |
---|
| 28 | |
---|
| 29 | #include "statistic_tests.hpp" |
---|
| 30 | #include "integrate.hpp" |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | namespace boost { |
---|
| 34 | namespace random { |
---|
| 35 | |
---|
| 36 | // Wikramaratna 1989 ACORN |
---|
| 37 | template<class IntType, int k, IntType m, IntType val> |
---|
| 38 | class additive_congruential |
---|
| 39 | { |
---|
| 40 | public: |
---|
| 41 | typedef IntType result_type; |
---|
| 42 | #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION |
---|
| 43 | static const bool has_fixed_range = true; |
---|
| 44 | static const result_type min_value = 0; |
---|
| 45 | static const result_type max_value = m-1; |
---|
| 46 | #else |
---|
| 47 | enum { |
---|
| 48 | has_fixed_range = true, |
---|
| 49 | min_value = 0, |
---|
| 50 | max_value = m-1 |
---|
| 51 | }; |
---|
| 52 | #endif |
---|
| 53 | template<class InputIterator> |
---|
| 54 | explicit additive_congruential(InputIterator start) { seed(start); } |
---|
| 55 | template<class InputIterator> |
---|
| 56 | void seed(InputIterator start) |
---|
| 57 | { |
---|
| 58 | for(int i = 0; i <= k; ++i, ++start) |
---|
| 59 | values[i] = *start; |
---|
| 60 | } |
---|
| 61 | |
---|
| 62 | result_type operator()() |
---|
| 63 | { |
---|
| 64 | for(int i = 1; i <= k; ++i) { |
---|
| 65 | IntType tmp = values[i-1] + values[i]; |
---|
| 66 | if(tmp >= m) |
---|
| 67 | tmp -= m; |
---|
| 68 | values[i] = tmp; |
---|
| 69 | } |
---|
| 70 | return values[k]; |
---|
| 71 | } |
---|
| 72 | result_type validation() const { return val; } |
---|
| 73 | private: |
---|
| 74 | IntType values[k+1]; |
---|
| 75 | }; |
---|
| 76 | |
---|
| 77 | |
---|
| 78 | template<class IntType, int r, int s, IntType m, IntType val> |
---|
| 79 | class lagged_fibonacci_int |
---|
| 80 | { |
---|
| 81 | public: |
---|
| 82 | typedef IntType result_type; |
---|
| 83 | #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION |
---|
| 84 | static const bool has_fixed_range = true; |
---|
| 85 | static const result_type min_value = 0; |
---|
| 86 | static const result_type max_value = m-1; |
---|
| 87 | #else |
---|
| 88 | enum { |
---|
| 89 | has_fixed_range = true, |
---|
| 90 | min_value = 0, |
---|
| 91 | max_value = m-1 |
---|
| 92 | }; |
---|
| 93 | #endif |
---|
| 94 | explicit lagged_fibonacci_int(IntType start) { seed(start); } |
---|
| 95 | template<class Generator> |
---|
| 96 | explicit lagged_fibonacci_int(Generator & gen) { seed(gen); } |
---|
| 97 | void seed(IntType start) |
---|
| 98 | { |
---|
| 99 | linear_congruential<uint32_t, 299375077, 0, 0, 0> init; |
---|
| 100 | seed(init); |
---|
| 101 | } |
---|
| 102 | template<class Generator> |
---|
| 103 | void seed(Generator & gen) |
---|
| 104 | { |
---|
| 105 | assert(r > s); |
---|
| 106 | for(int i = 0; i < 607; ++i) |
---|
| 107 | values[i] = gen(); |
---|
| 108 | current = 0; |
---|
| 109 | lag = r-s; |
---|
| 110 | } |
---|
| 111 | |
---|
| 112 | result_type operator()() |
---|
| 113 | { |
---|
| 114 | result_type tmp = values[current] + values[lag]; |
---|
| 115 | if(tmp >= m) |
---|
| 116 | tmp -= m; |
---|
| 117 | values[current] = tmp; |
---|
| 118 | ++current; |
---|
| 119 | if(current >= r) |
---|
| 120 | current = 0; |
---|
| 121 | ++lag; |
---|
| 122 | if(lag >= r) |
---|
| 123 | lag = 0; |
---|
| 124 | return tmp; |
---|
| 125 | } |
---|
| 126 | result_type validation() const { return val; } |
---|
| 127 | private: |
---|
| 128 | result_type values[r]; |
---|
| 129 | int current, lag; |
---|
| 130 | }; |
---|
| 131 | |
---|
| 132 | } // namespace random |
---|
| 133 | } // namespace boost |
---|
| 134 | |
---|
| 135 | // distributions from Haertel's dissertation |
---|
| 136 | // (additional parameterizations of the basic templates) |
---|
| 137 | namespace Haertel { |
---|
| 138 | typedef boost::random::linear_congruential<boost::uint64_t, 45965, 453816691, |
---|
| 139 | (boost::uint64_t(1)<<31), 0> LCG_Af2; |
---|
| 140 | typedef boost::random::linear_congruential<boost::uint64_t, 211936855, 0, |
---|
| 141 | (boost::uint64_t(1)<<29)-3, 0> LCG_Die1; |
---|
| 142 | typedef boost::random::linear_congruential<boost::uint32_t, 2824527309u, 0, |
---|
| 143 | 0, 0> LCG_Fis; |
---|
| 144 | typedef boost::random::linear_congruential<boost::uint64_t, 950706376u, 0, |
---|
| 145 | (boost::uint64_t(1)<<31)-1, 0> LCG_FM; |
---|
| 146 | typedef boost::random::linear_congruential<boost::int32_t, 51081, 0, |
---|
| 147 | 2147483647, 0> LCG_Hae; |
---|
| 148 | typedef boost::random::linear_congruential<boost::uint32_t, 69069, 1, |
---|
| 149 | 0, 0> LCG_VAX; |
---|
| 150 | typedef boost::random::inversive_congruential<boost::int64_t, 240318, 197, |
---|
| 151 | 1000081, 0> NLG_Inv1; |
---|
| 152 | typedef boost::random::inversive_congruential<boost::int64_t, 15707262, |
---|
| 153 | 13262967, (1<<24)-17, 0> NLG_Inv2; |
---|
| 154 | typedef boost::random::inversive_congruential<boost::int32_t, 1, 1, |
---|
| 155 | 2147483647, 0> NLG_Inv4; |
---|
| 156 | typedef boost::random::inversive_congruential<boost::int32_t, 1, 2, |
---|
| 157 | 1<<30, 0> NLG_Inv5; |
---|
| 158 | typedef boost::random::additive_congruential<boost::int32_t, 6, |
---|
| 159 | (1<<30)-35, 0> MRG_Acorn7; |
---|
| 160 | typedef boost::random::lagged_fibonacci_int<boost::uint32_t, 607, 273, |
---|
| 161 | 0, 0> MRG_Fib2; |
---|
| 162 | |
---|
| 163 | template<class Gen, class T> |
---|
| 164 | inline void check_validation(Gen & gen, T value, const std::string & name) |
---|
| 165 | { |
---|
| 166 | for(int i = 0; i < 100000-1; ++i) |
---|
| 167 | gen(); |
---|
| 168 | if(value != gen()) |
---|
| 169 | std::cout << name << ": validation failed" << std::endl; |
---|
| 170 | } |
---|
| 171 | |
---|
| 172 | // we have validation after 100000 steps with Haertel's generators |
---|
| 173 | template<class Gen, class T> |
---|
| 174 | void validate(T value, const std::string & name) |
---|
| 175 | { |
---|
| 176 | Gen gen(1234567); |
---|
| 177 | check_validation(gen, value, name); |
---|
| 178 | } |
---|
| 179 | |
---|
| 180 | void validate_all() |
---|
| 181 | { |
---|
| 182 | validate<LCG_Af2>(183269031u, "LCG_Af2"); |
---|
| 183 | validate<LCG_Die1>(522319944u, "LCG_Die1"); |
---|
| 184 | validate<LCG_Fis>(-2065162233u, "LCG_Fis"); |
---|
| 185 | validate<LCG_FM>(581815473u, "LCG_FM"); |
---|
| 186 | validate<LCG_Hae>(28931709, "LCG_Hae"); |
---|
| 187 | validate<LCG_VAX>(1508154087u, "LCG_VAX"); |
---|
| 188 | validate<NLG_Inv2>(6666884, "NLG_Inv2"); |
---|
| 189 | validate<NLG_Inv4>(1521640076, "NLG_Inv4"); |
---|
| 190 | validate<NLG_Inv5>(641840839, "NLG_Inv5"); |
---|
| 191 | static const int acorn7_init[] |
---|
| 192 | = { 1234567, 7654321, 246810, 108642, 13579, 97531, 555555 }; |
---|
| 193 | MRG_Acorn7 acorn7(acorn7_init); |
---|
| 194 | check_validation(acorn7, 874294697, "MRG_Acorn7"); |
---|
| 195 | validate<MRG_Fib2>(1234567u, "MRG_Fib2"); |
---|
| 196 | } |
---|
| 197 | } // namespace Haertel |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | |
---|
| 202 | double normal_density(double x) |
---|
| 203 | { |
---|
| 204 | const double pi = 3.14159265358979323846; |
---|
| 205 | return 1/std::sqrt(2*pi) * std::exp(-x*x/2); |
---|
| 206 | } |
---|
| 207 | |
---|
| 208 | namespace std { |
---|
| 209 | #ifdef _CXXRTCF_H__ |
---|
| 210 | using _CS_swamp::lgamma; |
---|
| 211 | #elif defined __SGI_STL_PORT |
---|
| 212 | using ::lgamma; |
---|
| 213 | #endif |
---|
| 214 | } |
---|
| 215 | |
---|
| 216 | |
---|
| 217 | class chi_square_density : public std::unary_function<double, double> |
---|
| 218 | { |
---|
| 219 | public: |
---|
| 220 | chi_square_density(int freedom) |
---|
| 221 | : _exponent( static_cast<double>(freedom)/2-1 ), |
---|
| 222 | _factor(1/(std::pow(2, _exponent+1) * std::exp(lgamma(_exponent+1)))) |
---|
| 223 | { } |
---|
| 224 | |
---|
| 225 | double operator()(double x) |
---|
| 226 | { |
---|
| 227 | return _factor*std::pow(x, _exponent)*std::exp(-x/2); |
---|
| 228 | } |
---|
| 229 | private: |
---|
| 230 | double _exponent, _factor; |
---|
| 231 | }; |
---|
| 232 | |
---|
| 233 | // computes F(x) or F(y) - F(x) |
---|
| 234 | class chi_square_probability : public distribution_function<double> |
---|
| 235 | { |
---|
| 236 | public: |
---|
| 237 | chi_square_probability(int freedom) : dens(freedom) {} |
---|
| 238 | double operator()(double x) { return operator()(0, x); } |
---|
| 239 | double operator()(double x, double y) |
---|
| 240 | { return trapezoid(dens, x, y, 1000); } |
---|
| 241 | private: |
---|
| 242 | chi_square_density dens; |
---|
| 243 | }; |
---|
| 244 | |
---|
| 245 | class uniform_distribution : public distribution_function<double> |
---|
| 246 | { |
---|
| 247 | public: |
---|
| 248 | uniform_distribution(double from, double to) : from(from), to(to) |
---|
| 249 | { assert(from < to); } |
---|
| 250 | double operator()(double x) |
---|
| 251 | { |
---|
| 252 | if(x < from) |
---|
| 253 | return 0; |
---|
| 254 | else if(x > to) |
---|
| 255 | return 1; |
---|
| 256 | else |
---|
| 257 | return (x-from)/(to-from); |
---|
| 258 | } |
---|
| 259 | double operator()(double x, double delta) |
---|
| 260 | { return operator()(x+delta) - operator()(x); } |
---|
| 261 | private: |
---|
| 262 | double from, to; |
---|
| 263 | }; |
---|
| 264 | |
---|
| 265 | class test_environment; |
---|
| 266 | |
---|
| 267 | class test_base |
---|
| 268 | { |
---|
| 269 | protected: |
---|
| 270 | explicit test_base(test_environment & env) : environment(env) { } |
---|
| 271 | void check(double val) const; |
---|
| 272 | private: |
---|
| 273 | test_environment & environment; |
---|
| 274 | }; |
---|
| 275 | |
---|
| 276 | class equidistribution_test : test_base |
---|
| 277 | { |
---|
| 278 | public: |
---|
| 279 | equidistribution_test(test_environment & env, unsigned int classes, |
---|
| 280 | unsigned int high_classes) |
---|
| 281 | : test_base(env), classes(classes), |
---|
| 282 | test_distrib_chi_square(chi_square_probability(classes-1), high_classes) |
---|
| 283 | { } |
---|
| 284 | |
---|
| 285 | template<class RNG> |
---|
| 286 | void run(RNG & rng, int n1, int n2) |
---|
| 287 | { |
---|
| 288 | using namespace boost; |
---|
| 289 | std::cout << "equidistribution: " << std::flush; |
---|
| 290 | equidistribution_experiment equi(classes); |
---|
| 291 | uniform_smallint<RNG> uint_linear(rng, 0, classes-1); |
---|
| 292 | check(run_experiment(test_distrib_chi_square, |
---|
| 293 | experiment_generator(equi, uint_linear, n1), n2)); |
---|
| 294 | check(run_experiment(test_distrib_chi_square, |
---|
| 295 | experiment_generator(equi, uint_linear, n1), 2*n2)); |
---|
| 296 | |
---|
| 297 | std::cout << " 2D: " << std::flush; |
---|
| 298 | equidistribution_2d_experiment equi_2d(classes); |
---|
| 299 | unsigned int root = static_cast<unsigned int>(std::sqrt(double(classes))); |
---|
| 300 | assert(root * root == classes); |
---|
| 301 | uniform_smallint<RNG> uint_square(rng, 0, root-1); |
---|
| 302 | check(run_experiment(test_distrib_chi_square, |
---|
| 303 | experiment_generator(equi_2d, uint_square, n1), n2)); |
---|
| 304 | check(run_experiment(test_distrib_chi_square, |
---|
| 305 | experiment_generator(equi_2d, uint_square, n1), 2*n2)); |
---|
| 306 | std::cout << std::endl; |
---|
| 307 | } |
---|
| 308 | private: |
---|
| 309 | unsigned int classes; |
---|
| 310 | distribution_experiment test_distrib_chi_square; |
---|
| 311 | }; |
---|
| 312 | |
---|
| 313 | class ks_equidistribution_test : test_base |
---|
| 314 | { |
---|
| 315 | public: |
---|
| 316 | ks_equidistribution_test(test_environment & env, unsigned int classes) |
---|
| 317 | : test_base(env), |
---|
| 318 | test_distrib_chi_square(kolmogorov_smirnov_probability(5000), |
---|
| 319 | classes) |
---|
| 320 | { } |
---|
| 321 | |
---|
| 322 | template<class RNG> |
---|
| 323 | void run(RNG & rng, int n1, int n2) |
---|
| 324 | { |
---|
| 325 | using namespace boost; |
---|
| 326 | std::cout << "KS: " << std::flush; |
---|
| 327 | // generator_reference_t<RNG> gen_ref(rng); |
---|
| 328 | RNG& gen_ref(rng); |
---|
| 329 | kolmogorov_experiment ks(n1); |
---|
| 330 | uniform_distribution ud((rng.min)(), (rng.max)()); |
---|
| 331 | check(run_experiment(test_distrib_chi_square, |
---|
| 332 | ks_experiment_generator(ks, gen_ref, ud), n2)); |
---|
| 333 | check(run_experiment(test_distrib_chi_square, |
---|
| 334 | ks_experiment_generator(ks, gen_ref, ud), 2*n2)); |
---|
| 335 | } |
---|
| 336 | private: |
---|
| 337 | distribution_experiment test_distrib_chi_square; |
---|
| 338 | }; |
---|
| 339 | |
---|
| 340 | class runs_test : test_base |
---|
| 341 | { |
---|
| 342 | public: |
---|
| 343 | runs_test(test_environment & env, unsigned int classes, |
---|
| 344 | unsigned int high_classes) |
---|
| 345 | : test_base(env), classes(classes), |
---|
| 346 | test_distrib_chi_square(chi_square_probability(classes-1), high_classes) |
---|
| 347 | { } |
---|
| 348 | |
---|
| 349 | template<class RNG> |
---|
| 350 | void run(RNG & rng, int n1, int n2) |
---|
| 351 | { |
---|
| 352 | using namespace boost; |
---|
| 353 | std::cout << "runs: up: " << std::flush; |
---|
| 354 | runs_experiment<true> r_up(classes); |
---|
| 355 | // generator_reference_t<RNG> gen_ref(rng); |
---|
| 356 | RNG& gen_ref(rng); |
---|
| 357 | check(run_experiment(test_distrib_chi_square, |
---|
| 358 | experiment_generator(r_up, gen_ref, n1), n2)); |
---|
| 359 | check(run_experiment(test_distrib_chi_square, |
---|
| 360 | experiment_generator(r_up, gen_ref, n1), 2*n2)); |
---|
| 361 | |
---|
| 362 | std::cout << " down: " << std::flush; |
---|
| 363 | runs_experiment<false> r_down(classes); |
---|
| 364 | check(run_experiment(test_distrib_chi_square, |
---|
| 365 | experiment_generator(r_down, gen_ref, n1), n2)); |
---|
| 366 | check(run_experiment(test_distrib_chi_square, |
---|
| 367 | experiment_generator(r_down, gen_ref, n1), 2*n2)); |
---|
| 368 | |
---|
| 369 | std::cout << std::endl; |
---|
| 370 | } |
---|
| 371 | private: |
---|
| 372 | unsigned int classes; |
---|
| 373 | distribution_experiment test_distrib_chi_square; |
---|
| 374 | }; |
---|
| 375 | |
---|
| 376 | class gap_test : test_base |
---|
| 377 | { |
---|
| 378 | public: |
---|
| 379 | gap_test(test_environment & env, unsigned int classes, |
---|
| 380 | unsigned int high_classes) |
---|
| 381 | : test_base(env), classes(classes), |
---|
| 382 | test_distrib_chi_square(chi_square_probability(classes-1), high_classes) |
---|
| 383 | { } |
---|
| 384 | |
---|
| 385 | template<class RNG> |
---|
| 386 | void run(RNG & rng, int n1, int n2) |
---|
| 387 | { |
---|
| 388 | using namespace boost; |
---|
| 389 | std::cout << "gaps: " << std::flush; |
---|
| 390 | gap_experiment gap(classes, 0.2, 0.8); |
---|
| 391 | // generator_reference_t<RNG> gen_ref(rng); |
---|
| 392 | RNG& gen_ref(rng); |
---|
| 393 | check(run_experiment(test_distrib_chi_square, |
---|
| 394 | experiment_generator(gap, gen_ref, n1), n2)); |
---|
| 395 | check(run_experiment(test_distrib_chi_square, |
---|
| 396 | experiment_generator(gap, gen_ref, n1), 2*n2)); |
---|
| 397 | |
---|
| 398 | std::cout << std::endl; |
---|
| 399 | } |
---|
| 400 | private: |
---|
| 401 | unsigned int classes; |
---|
| 402 | distribution_experiment test_distrib_chi_square; |
---|
| 403 | }; |
---|
| 404 | |
---|
| 405 | class poker_test : test_base |
---|
| 406 | { |
---|
| 407 | public: |
---|
| 408 | poker_test(test_environment & env, unsigned int classes, |
---|
| 409 | unsigned int high_classes) |
---|
| 410 | : test_base(env), classes(classes), |
---|
| 411 | test_distrib_chi_square(chi_square_probability(classes-1), high_classes) |
---|
| 412 | { } |
---|
| 413 | |
---|
| 414 | template<class RNG> |
---|
| 415 | void run(RNG & rng, int n1, int n2) |
---|
| 416 | { |
---|
| 417 | using namespace boost; |
---|
| 418 | std::cout << "poker: " << std::flush; |
---|
| 419 | poker_experiment poker(8, classes); |
---|
| 420 | uniform_smallint<RNG> usmall(rng, 0, 7); |
---|
| 421 | check(run_experiment(test_distrib_chi_square, |
---|
| 422 | experiment_generator(poker, usmall, n1), n2)); |
---|
| 423 | check(run_experiment(test_distrib_chi_square, |
---|
| 424 | experiment_generator(poker, usmall, n1), 2*n2)); |
---|
| 425 | std::cout << std::endl; |
---|
| 426 | } |
---|
| 427 | private: |
---|
| 428 | unsigned int classes; |
---|
| 429 | distribution_experiment test_distrib_chi_square; |
---|
| 430 | }; |
---|
| 431 | |
---|
| 432 | class coupon_collector_test : test_base |
---|
| 433 | { |
---|
| 434 | public: |
---|
| 435 | coupon_collector_test(test_environment & env, unsigned int classes, |
---|
| 436 | unsigned int high_classes) |
---|
| 437 | : test_base(env), classes(classes), |
---|
| 438 | test_distrib_chi_square(chi_square_probability(classes-1), high_classes) |
---|
| 439 | { } |
---|
| 440 | |
---|
| 441 | template<class RNG> |
---|
| 442 | void run(RNG & rng, int n1, int n2) |
---|
| 443 | { |
---|
| 444 | using namespace boost; |
---|
| 445 | std::cout << "coupon collector: " << std::flush; |
---|
| 446 | coupon_collector_experiment coupon(5, classes); |
---|
| 447 | |
---|
| 448 | uniform_smallint<RNG> usmall(rng, 0, 4); |
---|
| 449 | check(run_experiment(test_distrib_chi_square, |
---|
| 450 | experiment_generator(coupon, usmall, n1), n2)); |
---|
| 451 | check(run_experiment(test_distrib_chi_square, |
---|
| 452 | experiment_generator(coupon, usmall, n1), 2*n2)); |
---|
| 453 | std::cout << std::endl; |
---|
| 454 | } |
---|
| 455 | private: |
---|
| 456 | unsigned int classes; |
---|
| 457 | distribution_experiment test_distrib_chi_square; |
---|
| 458 | }; |
---|
| 459 | |
---|
| 460 | class permutation_test : test_base |
---|
| 461 | { |
---|
| 462 | public: |
---|
| 463 | permutation_test(test_environment & env, unsigned int classes, |
---|
| 464 | unsigned int high_classes) |
---|
| 465 | : test_base(env), classes(classes), |
---|
| 466 | test_distrib_chi_square(chi_square_probability(fac<int>(classes)-1), |
---|
| 467 | high_classes) |
---|
| 468 | { } |
---|
| 469 | |
---|
| 470 | template<class RNG> |
---|
| 471 | void run(RNG & rng, int n1, int n2) |
---|
| 472 | { |
---|
| 473 | using namespace boost; |
---|
| 474 | std::cout << "permutation: " << std::flush; |
---|
| 475 | permutation_experiment perm(classes); |
---|
| 476 | |
---|
| 477 | // generator_reference_t<RNG> gen_ref(rng); |
---|
| 478 | RNG& gen_ref(rng); |
---|
| 479 | check(run_experiment(test_distrib_chi_square, |
---|
| 480 | experiment_generator(perm, gen_ref, n1), n2)); |
---|
| 481 | check(run_experiment(test_distrib_chi_square, |
---|
| 482 | experiment_generator(perm, gen_ref, n1), 2*n2)); |
---|
| 483 | std::cout << std::endl; |
---|
| 484 | } |
---|
| 485 | private: |
---|
| 486 | unsigned int classes; |
---|
| 487 | distribution_experiment test_distrib_chi_square; |
---|
| 488 | }; |
---|
| 489 | |
---|
| 490 | class maximum_test : test_base |
---|
| 491 | { |
---|
| 492 | public: |
---|
| 493 | maximum_test(test_environment & env, unsigned int high_classes) |
---|
| 494 | : test_base(env), |
---|
| 495 | test_distrib_chi_square(kolmogorov_smirnov_probability(1000), |
---|
| 496 | high_classes) |
---|
| 497 | { } |
---|
| 498 | |
---|
| 499 | template<class RNG> |
---|
| 500 | void run(RNG & rng, int n1, int n2) |
---|
| 501 | { |
---|
| 502 | using namespace boost; |
---|
| 503 | std::cout << "maximum-of-t: " << std::flush; |
---|
| 504 | maximum_experiment<RNG> mx(rng, n1, 5); |
---|
| 505 | check(run_experiment(test_distrib_chi_square, mx, n2)); |
---|
| 506 | check(run_experiment(test_distrib_chi_square, mx, 2*n2)); |
---|
| 507 | std::cout << std::endl; |
---|
| 508 | } |
---|
| 509 | private: |
---|
| 510 | distribution_experiment test_distrib_chi_square; |
---|
| 511 | }; |
---|
| 512 | |
---|
| 513 | class birthday_test : test_base |
---|
| 514 | { |
---|
| 515 | public: |
---|
| 516 | birthday_test(test_environment & env) |
---|
| 517 | : test_base(env) |
---|
| 518 | { } |
---|
| 519 | |
---|
| 520 | template<class RNG> |
---|
| 521 | void run(RNG & rng, int n1, int n2) |
---|
| 522 | { |
---|
| 523 | using namespace boost; |
---|
| 524 | std::cout << "birthday spacing: " << std::flush; |
---|
| 525 | uniform_int<RNG> uni(rng, 0, (1<<25)-1); |
---|
| 526 | birthday_spacing_experiment bsp(4, 512, (1<<25)); |
---|
| 527 | std::cout << run_experiment(bsp, uni, n1); |
---|
| 528 | #if 0 |
---|
| 529 | check(run_experiment(test_distrib_chi_square, |
---|
| 530 | experiment_generator(perm, gen_ref, n1), n2)); |
---|
| 531 | check(run_experiment(test_distrib_chi_square, |
---|
| 532 | experiment_generator(perm, gen_ref, n1), 2*n2)); |
---|
| 533 | #endif |
---|
| 534 | std::cout << std::endl; |
---|
| 535 | } |
---|
| 536 | |
---|
| 537 | |
---|
| 538 | }; |
---|
| 539 | |
---|
| 540 | class test_environment |
---|
| 541 | { |
---|
| 542 | public: |
---|
| 543 | static const int classes = 20; |
---|
| 544 | explicit test_environment(double confid) |
---|
| 545 | : confidence(confid), |
---|
| 546 | confidence_chi_square_quantil(quantil(chi_square_density(classes-1), 0, confidence, 1e-4)), |
---|
| 547 | test_distrib_chi_square6(chi_square_probability(7-1), classes), |
---|
| 548 | ksequi_test(*this, classes), |
---|
| 549 | equi_test(*this, 100, classes), |
---|
| 550 | rns_test(*this, 7, classes), |
---|
| 551 | gp_test(*this, 7, classes), |
---|
| 552 | pk_test(*this, 5, classes), |
---|
| 553 | cpn_test(*this, 15, classes), |
---|
| 554 | perm_test(*this, 5, classes), |
---|
| 555 | max_test(*this, classes), |
---|
| 556 | bday_test(*this) |
---|
| 557 | { |
---|
| 558 | std::cout << "Confidence level: " << confid |
---|
| 559 | << "; 1-alpha = " << (1-confid) |
---|
| 560 | << "; chi_square(" << (classes-1) |
---|
| 561 | << ", " << confidence_chi_square_quantil |
---|
| 562 | << ") = " |
---|
| 563 | << chi_square_probability(classes-1)(0, confidence_chi_square_quantil) |
---|
| 564 | << std::endl; |
---|
| 565 | } |
---|
| 566 | |
---|
| 567 | bool check_confidence(double val, double chi_square_conf) const |
---|
| 568 | { |
---|
| 569 | std::cout << val; |
---|
| 570 | bool result = (val <= chi_square_conf); |
---|
| 571 | if(!result) { |
---|
| 572 | std::cout << "* ["; |
---|
| 573 | double prob = (val > 10*chi_square_conf ? 1 : |
---|
| 574 | chi_square_probability(classes-1)(0, val)); |
---|
| 575 | std::cout << (1-prob) << "]"; |
---|
| 576 | } |
---|
| 577 | std::cout << " " << std::flush; |
---|
| 578 | return result; |
---|
| 579 | } |
---|
| 580 | |
---|
| 581 | bool check(double chi_square_value) const |
---|
| 582 | { |
---|
| 583 | return check_confidence(chi_square_value, confidence_chi_square_quantil); |
---|
| 584 | } |
---|
| 585 | |
---|
| 586 | template<class RNG> |
---|
| 587 | void run_test(const std::string & name) |
---|
| 588 | { |
---|
| 589 | using namespace boost; |
---|
| 590 | |
---|
| 591 | std::cout << "Running tests on " << name << std::endl; |
---|
| 592 | |
---|
| 593 | RNG rng(1234567); |
---|
| 594 | typedef boost::uniform_01<RNG> UGen; |
---|
| 595 | |
---|
| 596 | #if 1 |
---|
| 597 | ksequi_test.run(rng, 5000, 250); |
---|
| 598 | equi_test.run(rng, 5000, 250); |
---|
| 599 | rns_test.run(rng, 100000, 250); |
---|
| 600 | gp_test.run(rng, 10000, 250); |
---|
| 601 | pk_test.run(rng, 5000, 250); |
---|
| 602 | cpn_test.run(rng, 500, 250); |
---|
| 603 | perm_test.run(rng, 1200, 250); |
---|
| 604 | max_test.run(rng, 1000, 250); |
---|
| 605 | #endif |
---|
| 606 | bday_test.run(rng, 1000, 150); |
---|
| 607 | |
---|
| 608 | std::cout << std::endl; |
---|
| 609 | } |
---|
| 610 | |
---|
| 611 | private: |
---|
| 612 | double confidence; |
---|
| 613 | double confidence_chi_square_quantil; |
---|
| 614 | distribution_experiment test_distrib_chi_square6; |
---|
| 615 | ks_equidistribution_test ksequi_test; |
---|
| 616 | equidistribution_test equi_test; |
---|
| 617 | runs_test rns_test; |
---|
| 618 | gap_test gp_test; |
---|
| 619 | poker_test pk_test; |
---|
| 620 | coupon_collector_test cpn_test; |
---|
| 621 | permutation_test perm_test; |
---|
| 622 | maximum_test max_test; |
---|
| 623 | birthday_test bday_test; |
---|
| 624 | }; |
---|
| 625 | |
---|
| 626 | void test_base::check(double val) const |
---|
| 627 | { |
---|
| 628 | environment.check(val); |
---|
| 629 | } |
---|
| 630 | |
---|
| 631 | void print_ks_table() |
---|
| 632 | { |
---|
| 633 | std::cout.setf(std::ios::fixed); |
---|
| 634 | std::cout.precision(5); |
---|
| 635 | static const double all_p[] = { 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99 }; |
---|
| 636 | for(int n = 0; n <= 10000; (n < 55 ? ++n : n *= 10)) { |
---|
| 637 | std::cout << std::setw(4) << n << " "; |
---|
| 638 | for(unsigned int i = 0; i < sizeof(all_p)/sizeof(all_p[0]); ++i) { |
---|
| 639 | std::cout << std::setw(8) |
---|
| 640 | << (n == 0 ? all_p[i] : |
---|
| 641 | invert_monotone_inc(kolmogorov_smirnov_probability(n), all_p[i], 0, 10)) |
---|
| 642 | << " "; |
---|
| 643 | } |
---|
| 644 | std::cout << std::endl; |
---|
| 645 | } |
---|
| 646 | } |
---|
| 647 | |
---|
| 648 | int main() |
---|
| 649 | { |
---|
| 650 | // Haertel::validate_all(); |
---|
| 651 | test_environment env(0.99); |
---|
| 652 | env.run_test<boost::minstd_rand>("minstd_rand"); |
---|
| 653 | env.run_test<boost::mt19937>("mt19937"); |
---|
| 654 | env.run_test<Haertel::LCG_Af2>("LCG_Af2"); |
---|
| 655 | env.run_test<Haertel::LCG_Die1>("LCG_Die1"); |
---|
| 656 | env.run_test<Haertel::LCG_Fis>("LCG_Fis"); |
---|
| 657 | env.run_test<Haertel::LCG_FM>("LCG_FM"); |
---|
| 658 | env.run_test<Haertel::LCG_Hae>("LCG_Hae"); |
---|
| 659 | env.run_test<Haertel::LCG_VAX>("LCG_VAX"); |
---|
| 660 | env.run_test<Haertel::NLG_Inv1>("NLG_Inv1"); |
---|
| 661 | env.run_test<Haertel::NLG_Inv2>("NLG_Inv2"); |
---|
| 662 | env.run_test<Haertel::NLG_Inv4>("NLG_Inv4"); |
---|
| 663 | env.run_test<Haertel::NLG_Inv5>("NLG_Inv5"); |
---|
| 664 | } |
---|