• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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
44bool 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>(max));
87    boost::variate_generator<boost::mt19937&, BOOST_RANDOM_DISTRIBUTION > vgen(gen, dist);
88
89    double prob = test.probability(test.run(vgen, 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
102template<class Dist1
103#ifdef BOOST_RANDOM_ARG2_NAME
104       , class Dist2
105#endif
106>
107bool 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
129int 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
142template<class T>
143bool handle_option(int& argc, char**& argv, const char* opt, T& value) {
144    if(std::strcmp(argv[0], 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
154int 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, "-r", repeat)
169             && !handle_option(argc, argv, "-" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG1_NAME), max_arg1)
170#ifdef BOOST_RANDOM_ARG2_TYPE
171             && !handle_option(argc, argv, "-" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG2_NAME), max_arg2)
172#endif
173             && !handle_option(argc, argv, "-t", 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