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/tools/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#include <sstream>
19#include <boost/archive/text_oarchive.hpp>
20#include <boost/archive/text_iarchive.hpp>
21
22using namespace boost;
23using namespace unit_test;
24using namespace boost::accumulators;
25
26///////////////////////////////////////////////////////////////////////////////
27// erf() not known by VC++ compiler!
28// my_erf() computes error function by numerically integrating with trapezoidal rule
29//
30double my_erf(double const& x, int const& n = 1000)
31{
32 double sum = 0.;
33 double delta = x/n;
34 for (int i = 1; i < n; ++i)
35 sum += std::exp(x: -i*i*delta*delta) * delta;
36 sum += 0.5 * delta * (1. + std::exp(x: -x*x));
37 return sum * 2. / std::sqrt(x: 3.141592653);
38}
39
40typedef accumulator_set<double, stats<tag::p_square_cumulative_distribution> > accumulator_t;
41typedef iterator_range<std::vector<std::pair<double, double> >::iterator > histogram_type;
42
43///////////////////////////////////////////////////////////////////////////////
44// test_stat
45//
46void test_stat()
47{
48 // tolerance in %
49 double epsilon = 3;
50
51 accumulator_t acc(p_square_cumulative_distribution_num_cells = 100);
52
53 // two random number generators
54 boost::lagged_fibonacci607 rng;
55 boost::normal_distribution<> mean_sigma(0,1);
56 boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal(rng, mean_sigma);
57
58 for (std::size_t i=0; i<1000000; ++i)
59 {
60 acc(normal());
61 }
62
63 histogram_type histogram = p_square_cumulative_distribution(acc);
64
65 for (std::size_t i = 0; i < histogram.size(); ++i)
66 {
67 // problem with small results: epsilon is relative (in percent), not absolute!
68 if ( histogram[i].second > 0.001 )
69 BOOST_CHECK_CLOSE( 0.5 * (1.0 + my_erf( histogram[i].first / std::sqrt(2.0) )), histogram[i].second, epsilon );
70 }
71}
72
73///////////////////////////////////////////////////////////////////////////////
74// test_persistency
75//
76void test_persistency()
77{
78 // "persistent" storage
79 std::stringstream ss;
80 // tolerance in %
81 double epsilon = 3;
82 {
83 accumulator_t acc(p_square_cumulative_distribution_num_cells = 100);
84
85 // two random number generators
86 boost::lagged_fibonacci607 rng;
87 boost::normal_distribution<> mean_sigma(0,1);
88 boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal(rng, mean_sigma);
89
90 for (std::size_t i=0; i<1000000; ++i)
91 {
92 acc(normal());
93 }
94
95 histogram_type histogram = p_square_cumulative_distribution(acc);
96
97 BOOST_CHECK_CLOSE(0.5 * (1.0 + my_erf(histogram[25].first / std::sqrt(2.0))), histogram[25].second, epsilon);
98 BOOST_CHECK_CLOSE(0.5 * (1.0 + my_erf(histogram[50].first / std::sqrt(2.0))), histogram[50].second, epsilon);
99 BOOST_CHECK_CLOSE(0.5 * (1.0 + my_erf(histogram[75].first / std::sqrt(2.0))), histogram[75].second, epsilon);
100 boost::archive::text_oarchive oa(ss);
101 acc.serialize(ar&: oa, file_version: 0);
102 }
103 accumulator_t acc(p_square_cumulative_distribution_num_cells = 100);
104 boost::archive::text_iarchive ia(ss);
105 acc.serialize(ar&: ia, file_version: 0);
106 histogram_type histogram = p_square_cumulative_distribution(acc);
107 BOOST_CHECK_CLOSE(0.5 * (1.0 + my_erf(histogram[25].first / std::sqrt(2.0))), histogram[25].second, epsilon);
108 BOOST_CHECK_CLOSE(0.5 * (1.0 + my_erf(histogram[50].first / std::sqrt(2.0))), histogram[50].second, epsilon);
109 BOOST_CHECK_CLOSE(0.5 * (1.0 + my_erf(histogram[75].first / std::sqrt(2.0))), histogram[75].second, epsilon);
110}
111
112///////////////////////////////////////////////////////////////////////////////
113// init_unit_test_suite
114//
115test_suite* init_unit_test_suite( int argc, char* argv[] )
116{
117 test_suite *test = BOOST_TEST_SUITE("p_square_cumulative_distribution test");
118
119 test->add(BOOST_TEST_CASE(&test_stat));
120 test->add(BOOST_TEST_CASE(&test_persistency));
121
122 return test;
123}
124
125

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