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/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/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 | // |
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(first: p1, last: p1 + sizeof(p1) / sizeof(double)); |
48 | |
49 | double p2[] = {0.001, 0.025}; |
50 | double p3[] = {0.975, 0.999}; |
51 | probs_normal1.assign(first: p2, last: p2 + sizeof(p2) / sizeof(double)); |
52 | probs_normal2.assign(first: p3, last: 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(first: p4, last: p4 + sizeof(p4) / sizeof(double)); |
57 | probs_normal_exact2.assign(first: p5, last: 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(x: -mu1 * (sample1 - 0.5 * mu1))); |
70 | acc_normal2(sample2, weight = std::exp(x: -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 | // |
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 | |