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