1// (C) Copyright Eric Niebler, Olivier Gygi 2006.
2// Use, modification and distribution are subject to the
3// Boost Software License, Version 1.0. (See accompanying file
4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6// Test case for p_square_cumul_dist.hpp
7
8#include <cmath>
9#include <boost/random.hpp>
10#include <boost/test/unit_test.hpp>
11#include <boost/test/floating_point_comparison.hpp>
12#include <boost/accumulators/numeric/functional/vector.hpp>
13#include <boost/accumulators/numeric/functional/complex.hpp>
14#include <boost/accumulators/numeric/functional/valarray.hpp>
15#include <boost/accumulators/accumulators.hpp>
16#include <boost/accumulators/statistics/stats.hpp>
17#include <boost/accumulators/statistics/p_square_cumul_dist.hpp>
18
19using namespace boost;
20using namespace unit_test;
21using namespace boost::accumulators;
22
23///////////////////////////////////////////////////////////////////////////////
24// erf() not known by VC++ compiler!
25// my_erf() computes error function by numerically integrating with trapezoidal rule
26//
27double my_erf(double const& x, int const& n = 1000)
28{
29 double sum = 0.;
30 double delta = x/n;
31 for (int i = 1; i < n; ++i)
32 sum += std::exp(x: -i*i*delta*delta) * delta;
33 sum += 0.5 * delta * (1. + std::exp(x: -x*x));
34 return sum * 2. / std::sqrt(x: 3.141592653);
35}
36
37///////////////////////////////////////////////////////////////////////////////
38// test_stat
39//
40void test_stat()
41{
42 // tolerance in %
43 double epsilon = 3;
44
45 typedef accumulator_set<double, stats<tag::p_square_cumulative_distribution> > accumulator_t;
46
47 accumulator_t acc(p_square_cumulative_distribution_num_cells = 100);
48
49 // two random number generators
50 boost::lagged_fibonacci607 rng;
51 boost::normal_distribution<> mean_sigma(0,1);
52 boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal(rng, mean_sigma);
53
54 for (std::size_t i=0; i<100000; ++i)
55 {
56 acc(normal());
57 }
58
59 typedef iterator_range<std::vector<std::pair<double, double> >::iterator > histogram_type;
60 histogram_type histogram = p_square_cumulative_distribution(acc);
61
62 for (std::size_t i = 0; i < histogram.size(); ++i)
63 {
64 // problem with small results: epsilon is relative (in percent), not absolute!
65 if ( histogram[i].second > 0.001 )
66 BOOST_CHECK_CLOSE( 0.5 * (1.0 + my_erf( histogram[i].first / std::sqrt(2.0) )), histogram[i].second, epsilon );
67 }
68}
69
70///////////////////////////////////////////////////////////////////////////////
71// init_unit_test_suite
72//
73test_suite* init_unit_test_suite( int argc, char* argv[] )
74{
75 test_suite *test = BOOST_TEST_SUITE("p_square_cumulative_distribution test");
76
77 test->add(BOOST_TEST_CASE(&test_stat));
78
79 return test;
80}
81
82

source code of boost/libs/accumulators/test/p_square_cumul_dist.cpp