1///////////////////////////////////////////////////////////////////////////////
2// extended_p_square_quantile.hpp
3//
4// Copyright 2005 Daniel Egloff. Distributed under the Boost
5// Software License, Version 1.0. (See accompanying file
6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
7
8#ifndef BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_QUANTILE_HPP_DE_01_01_2006
9#define BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_QUANTILE_HPP_DE_01_01_2006
10
11#include <vector>
12#include <functional>
13#include <boost/throw_exception.hpp>
14#include <boost/range/begin.hpp>
15#include <boost/range/end.hpp>
16#include <boost/range/iterator_range.hpp>
17#include <boost/iterator/transform_iterator.hpp>
18#include <boost/iterator/counting_iterator.hpp>
19#include <boost/iterator/permutation_iterator.hpp>
20#include <boost/parameter/keyword.hpp>
21#include <boost/mpl/placeholders.hpp>
22#include <boost/type_traits/is_same.hpp>
23#include <boost/accumulators/framework/accumulator_base.hpp>
24#include <boost/accumulators/framework/extractor.hpp>
25#include <boost/accumulators/numeric/functional.hpp>
26#include <boost/accumulators/framework/parameters/sample.hpp>
27#include <boost/accumulators/framework/depends_on.hpp>
28#include <boost/accumulators/statistics_fwd.hpp>
29#include <boost/accumulators/statistics/count.hpp>
30#include <boost/accumulators/statistics/parameters/quantile_probability.hpp>
31#include <boost/accumulators/statistics/extended_p_square.hpp>
32#include <boost/accumulators/statistics/weighted_extended_p_square.hpp>
33#include <boost/accumulators/statistics/times2_iterator.hpp>
34
35#ifdef _MSC_VER
36# pragma warning(push)
37# pragma warning(disable: 4127) // conditional expression is constant
38#endif
39
40namespace boost { namespace accumulators
41{
42
43namespace impl
44{
45 ///////////////////////////////////////////////////////////////////////////////
46 // extended_p_square_quantile_impl
47 // single quantile estimation
48 /**
49 @brief Quantile estimation using the extended \f$P^2\f$ algorithm for weighted and unweighted samples
50
51 Uses the quantile estimates calculated by the extended \f$P^2\f$ algorithm to compute
52 intermediate quantile estimates by means of quadratic interpolation.
53
54 @param quantile_probability The probability of the quantile to be estimated.
55 */
56 template<typename Sample, typename Impl1, typename Impl2> // Impl1: weighted/unweighted // Impl2: linear/quadratic
57 struct extended_p_square_quantile_impl
58 : accumulator_base
59 {
60 typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
61 typedef std::vector<float_type> array_type;
62 typedef iterator_range<
63 detail::lvalue_index_iterator<
64 permutation_iterator<
65 typename array_type::const_iterator
66 , detail::times2_iterator
67 >
68 >
69 > range_type;
70 // for boost::result_of
71 typedef float_type result_type;
72
73 template<typename Args>
74 extended_p_square_quantile_impl(Args const &args)
75 : probabilities(
76 boost::begin(args[extended_p_square_probabilities])
77 , boost::end(args[extended_p_square_probabilities])
78 )
79 , probability()
80 {
81 }
82
83 template<typename Args>
84 result_type result(Args const &args) const
85 {
86 typedef
87 typename mpl::if_<
88 is_same<Impl1, weighted>
89 , tag::weighted_extended_p_square
90 , tag::extended_p_square
91 >::type
92 extended_p_square_tag;
93
94 extractor<extended_p_square_tag> const some_extended_p_square = {};
95
96 array_type heights(some_extended_p_square(args).size());
97 std::copy(some_extended_p_square(args).begin(), some_extended_p_square(args).end(), heights.begin());
98
99 this->probability = args[quantile_probability];
100
101 typename array_type::const_iterator iter_probs = std::lower_bound(this->probabilities.begin(), this->probabilities.end(), this->probability);
102 std::size_t dist = std::distance(this->probabilities.begin(), iter_probs);
103 typename array_type::const_iterator iter_heights = heights.begin() + dist;
104
105 // If this->probability is not in a valid range return NaN or throw exception
106 if (this->probability < *this->probabilities.begin() || this->probability > *(this->probabilities.end() - 1))
107 {
108 if (std::numeric_limits<result_type>::has_quiet_NaN)
109 {
110 return std::numeric_limits<result_type>::quiet_NaN();
111 }
112 else
113 {
114 std::ostringstream msg;
115 msg << "probability = " << this->probability << " is not in valid range (";
116 msg << *this->probabilities.begin() << ", " << *(this->probabilities.end() - 1) << ")";
117 boost::throw_exception(e: std::runtime_error(msg.str()));
118 return Sample(0);
119 }
120
121 }
122
123 if (*iter_probs == this->probability)
124 {
125 return heights[dist];
126 }
127 else
128 {
129 result_type res;
130
131 if (is_same<Impl2, linear>::value)
132 {
133 /////////////////////////////////////////////////////////////////////////////////
134 // LINEAR INTERPOLATION
135 //
136 float_type p1 = *iter_probs;
137 float_type p0 = *(iter_probs - 1);
138 float_type h1 = *iter_heights;
139 float_type h0 = *(iter_heights - 1);
140
141 float_type a = numeric::fdiv(h1 - h0, p1 - p0);
142 float_type b = h1 - p1 * a;
143
144 res = a * this->probability + b;
145 }
146 else
147 {
148 /////////////////////////////////////////////////////////////////////////////////
149 // QUADRATIC INTERPOLATION
150 //
151 float_type p0, p1, p2;
152 float_type h0, h1, h2;
153
154 if ( (dist == 1 || *iter_probs - this->probability <= this->probability - *(iter_probs - 1) ) && dist != this->probabilities.size() - 1 )
155 {
156 p0 = *(iter_probs - 1);
157 p1 = *iter_probs;
158 p2 = *(iter_probs + 1);
159 h0 = *(iter_heights - 1);
160 h1 = *iter_heights;
161 h2 = *(iter_heights + 1);
162 }
163 else
164 {
165 p0 = *(iter_probs - 2);
166 p1 = *(iter_probs - 1);
167 p2 = *iter_probs;
168 h0 = *(iter_heights - 2);
169 h1 = *(iter_heights - 1);
170 h2 = *iter_heights;
171 }
172
173 float_type hp21 = numeric::fdiv(h2 - h1, p2 - p1);
174 float_type hp10 = numeric::fdiv(h1 - h0, p1 - p0);
175 float_type p21 = numeric::fdiv(p2 * p2 - p1 * p1, p2 - p1);
176 float_type p10 = numeric::fdiv(p1 * p1 - p0 * p0, p1 - p0);
177
178 float_type a = numeric::fdiv(hp21 - hp10, p21 - p10);
179 float_type b = hp21 - a * p21;
180 float_type c = h2 - a * p2 * p2 - b * p2;
181
182 res = a * this->probability * this-> probability + b * this->probability + c;
183 }
184
185 return res;
186 }
187
188 }
189
190 public:
191 // make this accumulator serializeable
192 // TODO: do we need to split to load/save and verify that the parameters did not change?
193 template<class Archive>
194 void serialize(Archive & ar, const unsigned int file_version)
195 {
196 ar & probabilities;
197 ar & probability;
198 }
199
200 private:
201
202 array_type probabilities;
203 mutable float_type probability;
204
205 };
206
207} // namespace impl
208
209///////////////////////////////////////////////////////////////////////////////
210// tag::extended_p_square_quantile
211//
212namespace tag
213{
214 struct extended_p_square_quantile
215 : depends_on<extended_p_square>
216 {
217 typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, unweighted, linear> impl;
218 };
219 struct extended_p_square_quantile_quadratic
220 : depends_on<extended_p_square>
221 {
222 typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, unweighted, quadratic> impl;
223 };
224 struct weighted_extended_p_square_quantile
225 : depends_on<weighted_extended_p_square>
226 {
227 typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, weighted, linear> impl;
228 };
229 struct weighted_extended_p_square_quantile_quadratic
230 : depends_on<weighted_extended_p_square>
231 {
232 typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, weighted, quadratic> impl;
233 };
234}
235
236///////////////////////////////////////////////////////////////////////////////
237// extract::extended_p_square_quantile
238// extract::weighted_extended_p_square_quantile
239//
240namespace extract
241{
242 extractor<tag::extended_p_square_quantile> const extended_p_square_quantile = {};
243 extractor<tag::extended_p_square_quantile_quadratic> const extended_p_square_quantile_quadratic = {};
244 extractor<tag::weighted_extended_p_square_quantile> const weighted_extended_p_square_quantile = {};
245 extractor<tag::weighted_extended_p_square_quantile_quadratic> const weighted_extended_p_square_quantile_quadratic = {};
246
247 BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_quantile)
248 BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_quantile_quadratic)
249 BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_extended_p_square_quantile)
250 BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_extended_p_square_quantile_quadratic)
251}
252
253using extract::extended_p_square_quantile;
254using extract::extended_p_square_quantile_quadratic;
255using extract::weighted_extended_p_square_quantile;
256using extract::weighted_extended_p_square_quantile_quadratic;
257
258// extended_p_square_quantile(linear) -> extended_p_square_quantile
259template<>
260struct as_feature<tag::extended_p_square_quantile(linear)>
261{
262 typedef tag::extended_p_square_quantile type;
263};
264
265// extended_p_square_quantile(quadratic) -> extended_p_square_quantile_quadratic
266template<>
267struct as_feature<tag::extended_p_square_quantile(quadratic)>
268{
269 typedef tag::extended_p_square_quantile_quadratic type;
270};
271
272// weighted_extended_p_square_quantile(linear) -> weighted_extended_p_square_quantile
273template<>
274struct as_feature<tag::weighted_extended_p_square_quantile(linear)>
275{
276 typedef tag::weighted_extended_p_square_quantile type;
277};
278
279// weighted_extended_p_square_quantile(quadratic) -> weighted_extended_p_square_quantile_quadratic
280template<>
281struct as_feature<tag::weighted_extended_p_square_quantile(quadratic)>
282{
283 typedef tag::weighted_extended_p_square_quantile_quadratic type;
284};
285
286// for the purposes of feature-based dependency resolution,
287// extended_p_square_quantile and weighted_extended_p_square_quantile
288// provide the same feature as quantile
289template<>
290struct feature_of<tag::extended_p_square_quantile>
291 : feature_of<tag::quantile>
292{
293};
294template<>
295struct feature_of<tag::extended_p_square_quantile_quadratic>
296 : feature_of<tag::quantile>
297{
298};
299// So that extended_p_square_quantile can be automatically substituted with
300// weighted_extended_p_square_quantile when the weight parameter is non-void
301template<>
302struct as_weighted_feature<tag::extended_p_square_quantile>
303{
304 typedef tag::weighted_extended_p_square_quantile type;
305};
306
307template<>
308struct feature_of<tag::weighted_extended_p_square_quantile>
309 : feature_of<tag::extended_p_square_quantile>
310{
311};
312
313// So that extended_p_square_quantile_quadratic can be automatically substituted with
314// weighted_extended_p_square_quantile_quadratic when the weight parameter is non-void
315template<>
316struct as_weighted_feature<tag::extended_p_square_quantile_quadratic>
317{
318 typedef tag::weighted_extended_p_square_quantile_quadratic type;
319};
320template<>
321struct feature_of<tag::weighted_extended_p_square_quantile_quadratic>
322 : feature_of<tag::extended_p_square_quantile_quadratic>
323{
324};
325
326}} // namespace boost::accumulators
327
328#ifdef _MSC_VER
329# pragma warning(pop)
330#endif
331
332#endif
333

source code of boost/libs/accumulators/include/boost/accumulators/statistics/extended_p_square_quantile.hpp