test_real_distribution.ipp

/* test_real_distribution.ipp
 *
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * $Id$
 *
 */

#ifndef BOOST_MATH_DISTRIBUTION_INIT
#ifdef BOOST_RANDOM_ARG2_TYPE
#define BOOST_MATH_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME, BOOST_RANDOM_ARG2_NAME)
#else
#define BOOST_MATH_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME)
#endif
#endif

#ifndef BOOST_RANDOM_DISTRIBUTION_INIT
#ifdef BOOST_RANDOM_ARG2_TYPE
#define BOOST_RANDOM_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME, BOOST_RANDOM_ARG2_NAME)
#else
#define BOOST_RANDOM_DISTRIBUTION_INIT (BOOST_RANDOM_ARG1_NAME)
#endif
#endif

#ifndef BOOST_RANDOM_P_CUTOFF
#define BOOST_RANDOM_P_CUTOFF 0.99
#endif

#include <boost/random/mersenne_twister.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/exception/diagnostic_information.hpp>
#include <boost/preprocessor/stringize.hpp>
#include <boost/range/numeric.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <iostream>
#include <vector>

#include "statistic_tests.hpp"
#include "chi_squared_test.hpp"

bool do_test(BOOST_RANDOM_ARG1_TYPE BOOST_RANDOM_ARG1_NAME,
#ifdef BOOST_RANDOM_ARG2_TYPE
             BOOST_RANDOM_ARG2_TYPE BOOST_RANDOM_ARG2_NAME,
#endif
             long long max, boost::mt19937& gen) {
    std::cout << "running " BOOST_PP_STRINGIZE(BOOST_RANDOM_DISTRIBUTION_NAME) "("
        << BOOST_RANDOM_ARG1_NAME;
#ifdef BOOST_RANDOM_ARG2_NAME
    std::cout << ", " << BOOST_RANDOM_ARG2_NAME;
#endif
    std::cout << ")" << " " << max << " times: " << std::flush;

    BOOST_MATH_DISTRIBUTION expected BOOST_MATH_DISTRIBUTION_INIT;
    
    BOOST_RANDOM_DISTRIBUTION dist BOOST_RANDOM_DISTRIBUTION_INIT;

#ifdef BOOST_RANDOM_DISTRIBUTION_MAX

    BOOST_RANDOM_DISTRIBUTION::result_type max_value = BOOST_RANDOM_DISTRIBUTION_MAX;

    std::vector<double> expected_pdf(max_value+1);
    {
        for(int i = 0; i <= max_value; ++i) {
            expected_pdf[i] = pdf(expected, i);
        }
        expected_pdf.back() += 1 - boost::accumulate(expected_pdf, 0.0);
    }
    
    std::vector<long long> results(max_value + 1);
    for(long long i = 0; i < max; ++i) {
        ++results[(std::min)(dist(gen), max_value)];
    }

    long long sum = boost::accumulate(results, 0ll);
    if(sum != max) {
        std::cout << "*** Failed: incorrect total: " << sum << " ***" << std::endl;
        return false;
    }
    double prob = chi_squared_test(results, expected_pdf, max);

#else

    kolmogorov_experiment test(boost::numeric_cast<int>(max));
    boost::variate_generator<boost::mt19937&, BOOST_RANDOM_DISTRIBUTION > vgen(gen, dist);

    double prob = test.probability(test.run(vgen, expected));

#endif

    bool result = prob < BOOST_RANDOM_P_CUTOFF;
    const char* err = result? "" : "*";
    std::cout << std::setprecision(17) << prob << err << std::endl;

    std::cout << std::setprecision(6);

    return result;
}

template<class Dist1
#ifdef BOOST_RANDOM_ARG2_NAME
       , class Dist2
#endif
>
bool do_tests(int repeat, Dist1 d1,
#ifdef BOOST_RANDOM_ARG2_NAME
              Dist2 d2,
#endif
              long long trials) {
    boost::mt19937 gen;
    int errors = 0;
    for(int i = 0; i < repeat; ++i) {
        if(!do_test(d1(gen),
#ifdef BOOST_RANDOM_ARG2_NAME
            d2(gen),
#endif
            trials, gen)) {
            ++errors;
        }
    }
    if(errors != 0) {
        std::cout << "*** " << errors << " errors detected ***" << std::endl;
    }
    return errors == 0;
}

int usage() {
    std::cerr << "Usage: test_" BOOST_PP_STRINGIZE(BOOST_RANDOM_DISTRIBUTION_NAME)
        " -r <repeat>"
        " -" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG1_NAME)
        " <max " BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG1_NAME) ">"
#ifdef BOOST_RANDOM_ARG2_NAME
        " -" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG2_NAME)
        " <max " BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG2_NAME) ">"
#endif
        " -t <trials>" << std::endl;
    return 2;
}

template<class T>
bool handle_option(int& argc, char**& argv, const char* opt, T& value) {
    if(std::strcmp(argv[0], opt) == 0 && argc > 1) {
        --argc;
        ++argv;
        value = boost::lexical_cast<T>(argv[0]);
        return true;
    } else {
        return false;
    }
}

int main(int argc, char** argv) {
    int repeat = 1;
    BOOST_RANDOM_ARG1_TYPE max_arg1 = BOOST_RANDOM_ARG1_DEFAULT;
#ifdef BOOST_RANDOM_ARG2_TYPE
    BOOST_RANDOM_ARG2_TYPE max_arg2 = BOOST_RANDOM_ARG2_DEFAULT;
#endif
    long long trials = 100000;

    if(argc > 0) {
        --argc;
        ++argv;
    }
    while(argc > 0) {
        if(argv[0][0] != '-') return usage();
        else if(!handle_option(argc, argv, "-r", repeat)
             && !handle_option(argc, argv, "-" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG1_NAME), max_arg1)
#ifdef BOOST_RANDOM_ARG2_TYPE
             && !handle_option(argc, argv, "-" BOOST_PP_STRINGIZE(BOOST_RANDOM_ARG2_NAME), max_arg2)
#endif
             && !handle_option(argc, argv, "-t", trials)) {
            return usage();
        }
        --argc;
        ++argv;
    }

    try {
        if(do_tests(repeat,
                    BOOST_RANDOM_ARG1_DISTRIBUTION(max_arg1),
#ifdef BOOST_RANDOM_ARG2_TYPE
                    BOOST_RANDOM_ARG2_DISTRIBUTION(max_arg2),
#endif
                    trials)) {
            return 0;
        } else {
            return EXIT_FAILURE;
        }
    } catch(...) {
        std::cerr << boost::current_exception_diagnostic_information() << std::endl;
        return EXIT_FAILURE;
    }
}