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_p_square_quantile.hpp |
7 | |
8 | #include <cmath> // for std::exp() |
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_p_square_quantile.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_p_square_quantile>, double> accumulator_t; |
29 | |
30 | // tolerance in % |
31 | double epsilon = 1; |
32 | |
33 | // some random number generators |
34 | double mu4 = -1.0; |
35 | double mu5 = -1.0; |
36 | double mu6 = 1.0; |
37 | double mu7 = 1.0; |
38 | boost::lagged_fibonacci607 rng; |
39 | boost::normal_distribution<> mean_sigma4(mu4, 1); |
40 | boost::normal_distribution<> mean_sigma5(mu5, 1); |
41 | boost::normal_distribution<> mean_sigma6(mu6, 1); |
42 | boost::normal_distribution<> mean_sigma7(mu7, 1); |
43 | boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal4(rng, mean_sigma4); |
44 | boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal5(rng, mean_sigma5); |
45 | boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal6(rng, mean_sigma6); |
46 | boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > normal7(rng, mean_sigma7); |
47 | |
48 | accumulator_t acc0(quantile_probability = 0.001); |
49 | accumulator_t acc1(quantile_probability = 0.025); |
50 | accumulator_t acc2(quantile_probability = 0.975); |
51 | accumulator_t acc3(quantile_probability = 0.999); |
52 | |
53 | accumulator_t acc4(quantile_probability = 0.001); |
54 | accumulator_t acc5(quantile_probability = 0.025); |
55 | accumulator_t acc6(quantile_probability = 0.975); |
56 | accumulator_t acc7(quantile_probability = 0.999); |
57 | |
58 | |
59 | for (std::size_t i=0; i<100000; ++i) |
60 | { |
61 | double sample = rng(); |
62 | acc0(sample, weight = 1.); |
63 | acc1(sample, weight = 1.); |
64 | acc2(sample, weight = 1.); |
65 | acc3(sample, weight = 1.); |
66 | |
67 | double sample4 = normal4(); |
68 | double sample5 = normal5(); |
69 | double sample6 = normal6(); |
70 | double sample7 = normal7(); |
71 | acc4(sample4, weight = std::exp(x: -mu4 * (sample4 - 0.5 * mu4))); |
72 | acc5(sample5, weight = std::exp(x: -mu5 * (sample5 - 0.5 * mu5))); |
73 | acc6(sample6, weight = std::exp(x: -mu6 * (sample6 - 0.5 * mu6))); |
74 | acc7(sample7, weight = std::exp(x: -mu7 * (sample7 - 0.5 * mu7))); |
75 | } |
76 | // check for uniform distribution with weight = 1 |
77 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc0), 0.001, 28 ); |
78 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc1), 0.025, 5 ); |
79 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc2), 0.975, epsilon ); |
80 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc3), 0.999, epsilon ); |
81 | |
82 | // check for shifted standard normal distribution ("importance sampling") |
83 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc4), -3.090232, epsilon ); |
84 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc5), -1.959963, epsilon ); |
85 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc6), 1.959963, epsilon ); |
86 | BOOST_CHECK_CLOSE( weighted_p_square_quantile(acc7), 3.090232, epsilon ); |
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_p_square_quantile test" ); |
95 | |
96 | test->add(BOOST_TEST_CASE(&test_stat)); |
97 | |
98 | return test; |
99 | } |
100 | |
101 | |