1 | /* test_real_distribution.ipp |
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 | #ifndef BOOST_MATH_DISTRIBUTION_INIT |
13 | #ifdef BOOST_RANDOM_ARG2_TYPE |
14 | #define BOOST_MATH_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME, BOOST_RANDOM_ARG2_NAME) |
15 | #else |
16 | #define BOOST_MATH_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME) |
17 | #endif |
18 | #endif |
19 | |
20 | #ifndef BOOST_RANDOM_DISTRIBUTION_INIT |
21 | #ifdef BOOST_RANDOM_ARG2_TYPE |
22 | #define BOOST_RANDOM_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME, BOOST_RANDOM_ARG2_NAME) |
23 | #else |
24 | #define BOOST_RANDOM_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME) |
25 | #endif |
26 | #endif |
27 | |
28 | #ifndef BOOST_RANDOM_P_CUTOFF |
29 | #define BOOST_RANDOM_P_CUTOFF 0.99 |
30 | #endif |
31 | |
32 | #include <boost/random/mersenne_twister.hpp> |
33 | #include <boost/lexical_cast.hpp> |
34 | #include <boost/exception/diagnostic_information.hpp> |
35 | #include <boost/preprocessor/stringize.hpp> |
36 | #include <boost/range/numeric.hpp> |
37 | #include <boost/numeric/conversion/cast.hpp> |
38 | #include <iostream> |
39 | #include <vector> |
40 | |
41 | #include "statistic_tests.hpp" |
42 | #include "chi_squared_test.hpp" |
43 | |
44 | bool do_test(BOOST_RANDOM_ARG1_TYPE BOOST_RANDOM_ARG1_NAME, |
45 | #ifdef BOOST_RANDOM_ARG2_TYPE |
46 | BOOST_RANDOM_ARG2_TYPE BOOST_RANDOM_ARG2_NAME, |
47 | #endif |
48 | long long max, boost::mt19937& gen) { |
49 | std::cout << "running " BOOST_PP_STRINGIZE(BOOST_RANDOM_DISTRIBUTION_NAME) "(" |
50 | << BOOST_RANDOM_ARG1_NAME; |
51 | #ifdef BOOST_RANDOM_ARG2_NAME |
52 | std::cout << ", " << BOOST_RANDOM_ARG2_NAME; |
53 | #endif |
54 | std::cout << ")" << " " << max << " times: " << std::flush; |
55 | |
56 | BOOST_MATH_DISTRIBUTION expected BOOST_MATH_DISTRIBUTION_INIT; |
57 | |
58 | BOOST_RANDOM_DISTRIBUTION dist BOOST_RANDOM_DISTRIBUTION_INIT; |
59 | |
60 | #ifdef BOOST_RANDOM_DISTRIBUTION_MAX |
61 | |
62 | BOOST_RANDOM_DISTRIBUTION::result_type max_value = BOOST_RANDOM_DISTRIBUTION_MAX; |
63 | |
64 | std::vector<double> expected_pdf(max_value+1); |
65 | { |
66 | for(int i = 0; i <= max_value; ++i) { |
67 | expected_pdf[i] = pdf(expected, i); |
68 | } |
69 | expected_pdf.back() += 1 - boost::accumulate(expected_pdf, 0.0); |
70 | } |
71 | |
72 | std::vector<long long> results(max_value + 1); |
73 | for(long long i = 0; i < max; ++i) { |
74 | ++results[(std::min)(dist(gen), max_value)]; |
75 | } |
76 | |
77 | long long sum = boost::accumulate(results, 0ll); |
78 | if(sum != max) { |
79 | std::cout << "*** Failed: incorrect total: " << sum << " ***" << std::endl; |
80 | return false; |
81 | } |
82 | double prob = chi_squared_test(results, expected_pdf, max); |
83 | |
84 | #else |
85 | |
86 | kolmogorov_experiment test(boost::numeric_cast<int>(arg: max)); |
87 | boost::variate_generator<boost::mt19937&, BOOST_RANDOM_DISTRIBUTION > vgen(gen, dist); |
88 | |
89 | double prob = test.probability(x: test.run(gen&: vgen, distrib: expected)); |
90 | |
91 | #endif |
92 | |
93 | bool result = prob < BOOST_RANDOM_P_CUTOFF; |
94 | const char* err = result? "" : "*" ; |
95 | std::cout << std::setprecision(17) << prob << err << std::endl; |
96 | |
97 | std::cout << std::setprecision(6); |
98 | |
99 | return result; |
100 | } |
101 | |
102 | template<class Dist1 |
103 | #ifdef BOOST_RANDOM_ARG2_NAME |
104 | , class Dist2 |
105 | #endif |
106 | > |
107 | bool do_tests(int repeat, Dist1 d1, |
108 | #ifdef BOOST_RANDOM_ARG2_NAME |
109 | Dist2 d2, |
110 | #endif |
111 | long long trials) { |
112 | boost::mt19937 gen; |
113 | int errors = 0; |
114 | for(int i = 0; i < repeat; ++i) { |
115 | if(!do_test(d1(gen), |
116 | #ifdef BOOST_RANDOM_ARG2_NAME |
117 | d2(gen), |
118 | #endif |
119 | trials, gen)) { |
120 | ++errors; |
121 | } |
122 | } |
123 | if(errors != 0) { |
124 | std::cout << "*** " << errors << " errors detected ***" << std::endl; |
125 | } |
126 | return errors == 0; |
127 | } |
128 | |
129 | int usage() { |
130 | std::cerr << "Usage: test_" BOOST_PP_STRINGIZE(BOOST_RANDOM_DISTRIBUTION_NAME) |
131 | " -r <repeat>" |
132 | " -" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG1_NAME) |
133 | " <max " BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG1_NAME) ">" |
134 | #ifdef BOOST_RANDOM_ARG2_NAME |
135 | " -" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG2_NAME) |
136 | " <max " BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG2_NAME) ">" |
137 | #endif |
138 | " -t <trials>" << std::endl; |
139 | return 2; |
140 | } |
141 | |
142 | template<class T> |
143 | bool handle_option(int& argc, char**& argv, const char* opt, T& value) { |
144 | if(std::strcmp(s1: argv[0], s2: opt) == 0 && argc > 1) { |
145 | --argc; |
146 | ++argv; |
147 | value = boost::lexical_cast<T>(argv[0]); |
148 | return true; |
149 | } else { |
150 | return false; |
151 | } |
152 | } |
153 | |
154 | int main(int argc, char** argv) { |
155 | int repeat = 1; |
156 | BOOST_RANDOM_ARG1_TYPE max_arg1 = BOOST_RANDOM_ARG1_DEFAULT; |
157 | #ifdef BOOST_RANDOM_ARG2_TYPE |
158 | BOOST_RANDOM_ARG2_TYPE max_arg2 = BOOST_RANDOM_ARG2_DEFAULT; |
159 | #endif |
160 | long long trials = 100000; |
161 | |
162 | if(argc > 0) { |
163 | --argc; |
164 | ++argv; |
165 | } |
166 | while(argc > 0) { |
167 | if(argv[0][0] != '-') return usage(); |
168 | else if(!handle_option(argc, argv, opt: "-r" , value&: repeat) |
169 | && !handle_option(argc, argv, opt: "-" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG1_NAME), value&: max_arg1) |
170 | #ifdef BOOST_RANDOM_ARG2_TYPE |
171 | && !handle_option(argc, argv, opt: "-" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG2_NAME), value&: max_arg2) |
172 | #endif |
173 | && !handle_option(argc, argv, opt: "-t" , value&: trials)) { |
174 | return usage(); |
175 | } |
176 | --argc; |
177 | ++argv; |
178 | } |
179 | |
180 | try { |
181 | if(do_tests(repeat, |
182 | BOOST_RANDOM_ARG1_DISTRIBUTION(max_arg1), |
183 | #ifdef BOOST_RANDOM_ARG2_TYPE |
184 | BOOST_RANDOM_ARG2_DISTRIBUTION(max_arg2), |
185 | #endif |
186 | trials)) { |
187 | return 0; |
188 | } else { |
189 | return EXIT_FAILURE; |
190 | } |
191 | } catch(...) { |
192 | std::cerr << boost::current_exception_diagnostic_information() << std::endl; |
193 | return EXIT_FAILURE; |
194 | } |
195 | } |
196 | |