1 // (C) Copyright Eric Niebler 2005.
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 weighted_extended_p_square.hpp
7
8 #include <iostream>
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/weighted_extended_p_square.hpp>
18
19 using namespace boost;
20 using namespace unit_test;
21 using namespace boost::accumulators;
22
23 ///////////////////////////////////////////////////////////////////////////////
24 // test_stat
25 //
test_stat()26 void test_stat()
27 {
28 typedef accumulator_set<double, stats<tag::weighted_extended_p_square>, double> accumulator_t;
29
30 // problem with small results: epsilon is relative (in percent), not absolute
31
32 // tolerance in %
33 double epsilon = 1;
34
35 // some random number generators
36 double mu1 = -1.0;
37 double mu2 = 1.0;
38 boost::lagged_fibonacci607 rng;
39 boost::normal_distribution<> mean_sigma1(mu1, 1);
40 boost::normal_distribution<> mean_sigma2(mu2, 1);
41 boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal1(rng, mean_sigma1);
42 boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal2(rng, mean_sigma2);
43
44 std::vector<double> probs_uniform, probs_normal1, probs_normal2, probs_normal_exact1, probs_normal_exact2;
45
46 double p1[] = {/*0.001,*/ 0.01, 0.1, 0.5, 0.9, 0.99, 0.999};
47 probs_uniform.assign(p1, p1 + sizeof(p1) / sizeof(double));
48
49 double p2[] = {0.001, 0.025};
50 double p3[] = {0.975, 0.999};
51 probs_normal1.assign(p2, p2 + sizeof(p2) / sizeof(double));
52 probs_normal2.assign(p3, p3 + sizeof(p3) / sizeof(double));
53
54 double p4[] = {-3.090232, -1.959963};
55 double p5[] = {1.959963, 3.090232};
56 probs_normal_exact1.assign(p4, p4 + sizeof(p4) / sizeof(double));
57 probs_normal_exact2.assign(p5, p5 + sizeof(p5) / sizeof(double));
58
59 accumulator_t acc_uniform(extended_p_square_probabilities = probs_uniform);
60 accumulator_t acc_normal1(extended_p_square_probabilities = probs_normal1);
61 accumulator_t acc_normal2(extended_p_square_probabilities = probs_normal2);
62
63 for (std::size_t i = 0; i < 100000; ++i)
64 {
65 acc_uniform(rng(), weight = 1.);
66
67 double sample1 = normal1();
68 double sample2 = normal2();
69 acc_normal1(sample1, weight = std::exp(-mu1 * (sample1 - 0.5 * mu1)));
70 acc_normal2(sample2, weight = std::exp(-mu2 * (sample2 - 0.5 * mu2)));
71 }
72
73 // check for uniform distribution
74 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_uniform)[0], probs_uniform[0], 6*epsilon);
75 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_uniform)[1], probs_uniform[1], 3*epsilon);
76 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_uniform)[2], probs_uniform[2], epsilon);
77 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_uniform)[3], probs_uniform[3], epsilon);
78 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_uniform)[4], probs_uniform[4], epsilon);
79 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_uniform)[5], probs_uniform[5], epsilon);
80
81 // check for standard normal distribution
82 for (std::size_t i = 0; i < probs_normal1.size(); ++i)
83 {
84 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_normal1)[i], probs_normal_exact1[i], epsilon);
85 BOOST_CHECK_CLOSE(weighted_extended_p_square(acc_normal2)[i], probs_normal_exact2[i], epsilon);
86 }
87 }
88
89 ///////////////////////////////////////////////////////////////////////////////
90 // init_unit_test_suite
91 //
init_unit_test_suite(int argc,char * argv[])92 test_suite* init_unit_test_suite( int argc, char* argv[] )
93 {
94 test_suite *test = BOOST_TEST_SUITE("weighted_extended_p_square test");
95
96 test->add(BOOST_TEST_CASE(&test_stat));
97
98 return test;
99 }
100
101