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 #include <sstream>
19 #include <boost/archive/text_oarchive.hpp>
20 #include <boost/archive/text_iarchive.hpp>
21
22 using namespace boost;
23 using namespace unit_test;
24 using 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 //
my_erf(double const & x,int const & n=1000)30 double 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(-i*i*delta*delta) * delta;
36 sum += 0.5 * delta * (1. + std::exp(-x*x));
37 return sum * 2. / std::sqrt(3.141592653);
38 }
39
40 typedef accumulator_set<double, stats<tag::p_square_cumulative_distribution> > accumulator_t;
41 typedef iterator_range<std::vector<std::pair<double, double> >::iterator > histogram_type;
42
43 ///////////////////////////////////////////////////////////////////////////////
44 // test_stat
45 //
test_stat()46 void 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 //
test_persistency()76 void 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(oa, 0);
102 }
103 accumulator_t acc(p_square_cumulative_distribution_num_cells = 100);
104 boost::archive::text_iarchive ia(ss);
105 acc.serialize(ia, 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 //
init_unit_test_suite(int argc,char * argv[])115 test_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