1 | /* test_piecewise_constant.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/piecewise_constant_distribution.hpp> |
13 | #include <boost/random/uniform_int.hpp> |
14 | #include <boost/random/mersenne_twister.hpp> |
15 | #include <boost/lexical_cast.hpp> |
16 | #include <boost/exception/diagnostic_information.hpp> |
17 | #include <boost/range/algorithm/lower_bound.hpp> |
18 | #include <boost/range/numeric.hpp> |
19 | #include <vector> |
20 | #include <iostream> |
21 | #include <iomanip> |
22 | |
23 | #include "statistic_tests.hpp" |
24 | |
25 | class piecewise_constant |
26 | { |
27 | public: |
28 | piecewise_constant(const std::vector<double>& intervals, const std::vector<double>& weights) |
29 | : intervals(intervals), |
30 | cumulative(1, 0.0) |
31 | { |
32 | boost::partial_sum(rng: weights, result: std::back_inserter(x&: cumulative)); |
33 | for(std::vector<double>::iterator iter = cumulative.begin(), end = cumulative.end(); |
34 | iter != end; ++iter) |
35 | { |
36 | *iter /= cumulative.back(); |
37 | } |
38 | } |
39 | |
40 | double cdf(double x) const |
41 | { |
42 | std::size_t index = boost::lower_bound(rng: intervals, val: x) - intervals.begin(); |
43 | if(index == 0) return 0; |
44 | else if(index == intervals.size()) return 1; |
45 | else { |
46 | double lower_weight = cumulative[index - 1]; |
47 | double upper_weight = cumulative[index]; |
48 | double lower = intervals[index - 1]; |
49 | double upper = intervals[index]; |
50 | return lower_weight + (x - lower) / (upper - lower) * (upper_weight - lower_weight); |
51 | } |
52 | } |
53 | private: |
54 | std::vector<double> intervals; |
55 | std::vector<double> cumulative; |
56 | }; |
57 | |
58 | double cdf(const piecewise_constant& dist, double x) |
59 | { |
60 | return dist.cdf(x); |
61 | } |
62 | |
63 | bool do_test(int n, int max) { |
64 | std::cout << "running piecewise_constant(p0, p1, ..., p" << n-1 << ")" << " " << max << " times: " << std::flush; |
65 | |
66 | std::vector<double> weights; |
67 | { |
68 | boost::mt19937 egen; |
69 | for(int i = 0; i < n; ++i) { |
70 | weights.push_back(x: egen()); |
71 | } |
72 | } |
73 | std::vector<double> intervals; |
74 | for(int i = 0; i <= n; ++i) { |
75 | intervals.push_back(x: i); |
76 | } |
77 | |
78 | piecewise_constant expected(intervals, weights); |
79 | |
80 | boost::random::piecewise_constant_distribution<> dist(intervals, weights); |
81 | boost::mt19937 gen; |
82 | kolmogorov_experiment test(max); |
83 | boost::variate_generator<boost::mt19937&, boost::random::piecewise_constant_distribution<> > vgen(gen, dist); |
84 | |
85 | double prob = test.probability(x: test.run(gen&: vgen, distrib: expected)); |
86 | |
87 | bool result = prob < 0.99; |
88 | const char* err = result? "" : "*" ; |
89 | std::cout << std::setprecision(17) << prob << err << std::endl; |
90 | |
91 | std::cout << std::setprecision(6); |
92 | |
93 | return result; |
94 | } |
95 | |
96 | bool do_tests(int repeat, int max_n, int trials) { |
97 | boost::mt19937 gen; |
98 | boost::uniform_int<> idist(1, max_n); |
99 | int errors = 0; |
100 | for(int i = 0; i < repeat; ++i) { |
101 | if(!do_test(n: idist(gen), max: trials)) { |
102 | ++errors; |
103 | } |
104 | } |
105 | if(errors != 0) { |
106 | std::cout << "*** " << errors << " errors detected ***" << std::endl; |
107 | } |
108 | return errors == 0; |
109 | } |
110 | |
111 | int usage() { |
112 | std::cerr << "Usage: test_piecewise_constant -r <repeat> -n <max n> -t <trials>" << std::endl; |
113 | return 2; |
114 | } |
115 | |
116 | template<class T> |
117 | bool handle_option(int& argc, char**& argv, char opt, T& value) { |
118 | if(argv[0][1] == opt && argc > 1) { |
119 | --argc; |
120 | ++argv; |
121 | value = boost::lexical_cast<T>(argv[0]); |
122 | return true; |
123 | } else { |
124 | return false; |
125 | } |
126 | } |
127 | |
128 | int main(int argc, char** argv) { |
129 | int repeat = 10; |
130 | int max_n = 10; |
131 | int trials = 1000000; |
132 | |
133 | if(argc > 0) { |
134 | --argc; |
135 | ++argv; |
136 | } |
137 | while(argc > 0) { |
138 | if(argv[0][0] != '-') return usage(); |
139 | else if(!handle_option(argc, argv, opt: 'r', value&: repeat) |
140 | && !handle_option(argc, argv, opt: 'n', value&: max_n) |
141 | && !handle_option(argc, argv, opt: 't', value&: trials)) { |
142 | return usage(); |
143 | } |
144 | --argc; |
145 | ++argv; |
146 | } |
147 | |
148 | try { |
149 | if(do_tests(repeat, max_n, trials)) { |
150 | return 0; |
151 | } else { |
152 | return EXIT_FAILURE; |
153 | } |
154 | } catch(...) { |
155 | std::cerr << boost::current_exception_diagnostic_information() << std::endl; |
156 | return EXIT_FAILURE; |
157 | } |
158 | } |
159 | |