| 1 | /* test_bernoulli.cpp |
| 2 | * |
| 3 | * Copyright Steven Watanabe 2011 |
| 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 | */ |
| 11 | |
| 12 | #include <boost/random/bernoulli_distribution.hpp> |
| 13 | #include <boost/random/uniform_int.hpp> |
| 14 | #include <boost/random/uniform_01.hpp> |
| 15 | #include <boost/random/mersenne_twister.hpp> |
| 16 | #include <boost/math/distributions/binomial.hpp> |
| 17 | #include <boost/lexical_cast.hpp> |
| 18 | #include <boost/exception/diagnostic_information.hpp> |
| 19 | #include <vector> |
| 20 | #include <iostream> |
| 21 | #include <numeric> |
| 22 | |
| 23 | #include "chi_squared_test.hpp" |
| 24 | |
| 25 | bool do_test(double p, long long max) { |
| 26 | std::cout << "running bernoulli(" << p << ")" << " " << max << " times: " << std::flush; |
| 27 | |
| 28 | boost::math::binomial expected(static_cast<double>(max), p); |
| 29 | |
| 30 | boost::random::bernoulli_distribution<> dist(p); |
| 31 | boost::mt19937 gen; |
| 32 | long long count = 0; |
| 33 | for(long long i = 0; i < max; ++i) { |
| 34 | if(dist(gen)) ++count; |
| 35 | } |
| 36 | |
| 37 | double prob = cdf(dist: expected, x: count); |
| 38 | |
| 39 | bool result = prob < 0.99 && prob > 0.01; |
| 40 | const char* err = result? "" : "*" ; |
| 41 | std::cout << std::setprecision(17) << prob << err << std::endl; |
| 42 | |
| 43 | std::cout << std::setprecision(6); |
| 44 | |
| 45 | return result; |
| 46 | } |
| 47 | |
| 48 | bool do_tests(int repeat, long long trials) { |
| 49 | boost::mt19937 gen; |
| 50 | boost::uniform_01<> rdist; |
| 51 | int errors = 0; |
| 52 | for(int i = 0; i < repeat; ++i) { |
| 53 | if(!do_test(p: rdist(gen), max: trials)) { |
| 54 | ++errors; |
| 55 | } |
| 56 | } |
| 57 | if(errors != 0) { |
| 58 | std::cout << "*** " << errors << " errors detected ***" << std::endl; |
| 59 | } |
| 60 | return errors == 0; |
| 61 | } |
| 62 | |
| 63 | int usage() { |
| 64 | std::cerr << "Usage: test_bernoulli_distribution -r <repeat> -t <trials>" << std::endl; |
| 65 | return 2; |
| 66 | } |
| 67 | |
| 68 | template<class T> |
| 69 | bool handle_option(int& argc, char**& argv, char opt, T& value) { |
| 70 | if(argv[0][1] == opt && argc > 1) { |
| 71 | --argc; |
| 72 | ++argv; |
| 73 | value = boost::lexical_cast<T>(argv[0]); |
| 74 | return true; |
| 75 | } else { |
| 76 | return false; |
| 77 | } |
| 78 | } |
| 79 | |
| 80 | int main(int argc, char** argv) { |
| 81 | int repeat = 10; |
| 82 | long long trials = 1000000ll; |
| 83 | |
| 84 | if(argc > 0) { |
| 85 | --argc; |
| 86 | ++argv; |
| 87 | } |
| 88 | while(argc > 0) { |
| 89 | if(argv[0][0] != '-') return usage(); |
| 90 | else if(!handle_option(argc, argv, opt: 'r', value&: repeat) |
| 91 | && !handle_option(argc, argv, opt: 't', value&: trials)) { |
| 92 | return usage(); |
| 93 | } |
| 94 | --argc; |
| 95 | ++argv; |
| 96 | } |
| 97 | |
| 98 | try { |
| 99 | if(do_tests(repeat, trials)) { |
| 100 | return 0; |
| 101 | } else { |
| 102 | return EXIT_FAILURE; |
| 103 | } |
| 104 | } catch(...) { |
| 105 | std::cerr << boost::current_exception_diagnostic_information() << std::endl; |
| 106 | return EXIT_FAILURE; |
| 107 | } |
| 108 | } |
| 109 | |