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 | |
40 | namespace boost { namespace accumulators |
41 | { |
42 | |
43 | namespace 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 | // |
212 | namespace 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 | // |
240 | namespace extract |
241 | { |
242 | extractor<tag::extended_p_square_quantile> const = {}; |
243 | extractor<tag::extended_p_square_quantile_quadratic> const = {}; |
244 | extractor<tag::weighted_extended_p_square_quantile> const = {}; |
245 | extractor<tag::weighted_extended_p_square_quantile_quadratic> const = {}; |
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 | |
253 | using extract::extended_p_square_quantile; |
254 | using extract::extended_p_square_quantile_quadratic; |
255 | using extract::weighted_extended_p_square_quantile; |
256 | using extract::weighted_extended_p_square_quantile_quadratic; |
257 | |
258 | // extended_p_square_quantile(linear) -> extended_p_square_quantile |
259 | template<> |
260 | struct 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 |
266 | template<> |
267 | struct 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 |
273 | template<> |
274 | struct 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 |
280 | template<> |
281 | struct 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 |
289 | template<> |
290 | struct feature_of<tag::extended_p_square_quantile> |
291 | : feature_of<tag::quantile> |
292 | { |
293 | }; |
294 | template<> |
295 | struct 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 |
301 | template<> |
302 | struct as_weighted_feature<tag::extended_p_square_quantile> |
303 | { |
304 | typedef tag::weighted_extended_p_square_quantile type; |
305 | }; |
306 | |
307 | template<> |
308 | struct 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 |
315 | template<> |
316 | struct as_weighted_feature<tag::extended_p_square_quantile_quadratic> |
317 | { |
318 | typedef tag::weighted_extended_p_square_quantile_quadratic type; |
319 | }; |
320 | template<> |
321 | struct 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 | |