1 ///////////////////////////////////////////////////////////////////////////////
2 // extended_p_square.hpp
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)
8 #ifndef BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_HPP_DE_01_01_2006
9 #define BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_HPP_DE_01_01_2006
13 #include <boost/range/begin.hpp>
14 #include <boost/range/end.hpp>
15 #include <boost/range/iterator_range.hpp>
16 #include <boost/iterator/transform_iterator.hpp>
17 #include <boost/iterator/counting_iterator.hpp>
18 #include <boost/iterator/permutation_iterator.hpp>
19 #include <boost/parameter/keyword.hpp>
20 #include <boost/mpl/placeholders.hpp>
21 #include <boost/accumulators/accumulators_fwd.hpp>
22 #include <boost/accumulators/framework/extractor.hpp>
23 #include <boost/accumulators/numeric/functional.hpp>
24 #include <boost/accumulators/framework/parameters/sample.hpp>
25 #include <boost/accumulators/framework/depends_on.hpp>
26 #include <boost/accumulators/statistics_fwd.hpp>
27 #include <boost/accumulators/statistics/count.hpp>
28 #include <boost/accumulators/statistics/times2_iterator.hpp>
30 namespace boost { namespace accumulators
32 ///////////////////////////////////////////////////////////////////////////////
33 // probabilities named parameter
35 BOOST_PARAMETER_NESTED_KEYWORD(tag, extended_p_square_probabilities, probabilities)
37 BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_probabilities)
41 ///////////////////////////////////////////////////////////////////////////////
42 // extended_p_square_impl
43 // multiple quantile estimation
45 @brief Multiple quantile estimation with the extended \f$P^2\f$ algorithm
47 Extended \f$P^2\f$ algorithm for estimation of several quantiles without storing samples.
48 Assume that \f$m\f$ quantiles \f$\xi_{p_1}, \ldots, \xi_{p_m}\f$ are to be estimated.
49 Instead of storing the whole sample cumulative distribution, the algorithm maintains only
50 \f$m+2\f$ principal markers and \f$m+1\f$ middle markers, whose positions are updated
51 with each sample and whose heights are adjusted (if necessary) using a piecewise-parablic
52 formula. The heights of these central markers are the current estimates of the quantiles
53 and returned as an iterator range.
55 For further details, see
57 K. E. E. Raatikainen, Simultaneous estimation of several quantiles, Simulation, Volume 49,
58 Number 4 (October), 1986, p. 159-164.
60 The extended \f$ P^2 \f$ algorithm generalizes the \f$ P^2 \f$ algorithm of
62 R. Jain and I. Chlamtac, The P^2 algorithm for dynamic calculation of quantiles and
63 histograms without storing observations, Communications of the ACM,
64 Volume 28 (October), Number 10, 1985, p. 1076-1085.
66 @param extended_p_square_probabilities A vector of quantile probabilities.
68 template<typename Sample>
69 struct extended_p_square_impl
72 typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
73 typedef std::vector<float_type> array_type;
74 // for boost::result_of
75 typedef iterator_range<
76 detail::lvalue_index_iterator<
78 typename array_type::const_iterator
79 , detail::times2_iterator
84 template<typename Args>
85 extended_p_square_impl(Args const &args)
87 boost::begin(args[extended_p_square_probabilities])
88 , boost::end(args[extended_p_square_probabilities])
90 , heights(2 * probabilities.size() + 3)
91 , actual_positions(heights.size())
92 , desired_positions(heights.size())
93 , positions_increments(heights.size())
95 std::size_t num_quantiles = this->probabilities.size();
96 std::size_t num_markers = this->heights.size();
98 for(std::size_t i = 0; i < num_markers; ++i)
100 this->actual_positions[i] = i + 1;
103 this->positions_increments[0] = 0.;
104 this->positions_increments[num_markers - 1] = 1.;
106 for(std::size_t i = 0; i < num_quantiles; ++i)
108 this->positions_increments[2 * i + 2] = probabilities[i];
111 for(std::size_t i = 0; i <= num_quantiles; ++i)
113 this->positions_increments[2 * i + 1] =
114 0.5 * (this->positions_increments[2 * i] + this->positions_increments[2 * i + 2]);
117 for(std::size_t i = 0; i < num_markers; ++i)
119 this->desired_positions[i] = 1. + 2. * (num_quantiles + 1.) * this->positions_increments[i];
123 template<typename Args>
124 void operator ()(Args const &args)
126 std::size_t cnt = count(args);
128 // m+2 principal markers and m+1 middle markers
129 std::size_t num_markers = 2 * this->probabilities.size() + 3;
131 // first accumulate num_markers samples
132 if(cnt <= num_markers)
134 this->heights[cnt - 1] = args[sample];
136 // complete the initialization of heights by sorting
137 if(cnt == num_markers)
139 std::sort(this->heights.begin(), this->heights.end());
144 std::size_t sample_cell = 1;
146 // find cell k = sample_cell such that heights[k-1] <= sample < heights[k]
147 if(args[sample] < this->heights[0])
149 this->heights[0] = args[sample];
152 else if(args[sample] >= this->heights[num_markers - 1])
154 this->heights[num_markers - 1] = args[sample];
155 sample_cell = num_markers - 1;
159 typedef typename array_type::iterator iterator;
160 iterator it = std::upper_bound(
161 this->heights.begin()
162 , this->heights.end()
166 sample_cell = std::distance(this->heights.begin(), it);
169 // update actual positions of all markers above sample_cell index
170 for(std::size_t i = sample_cell; i < num_markers; ++i)
172 ++this->actual_positions[i];
175 // update desired positions of all markers
176 for(std::size_t i = 0; i < num_markers; ++i)
178 this->desired_positions[i] += this->positions_increments[i];
181 // adjust heights and actual positions of markers 1 to num_markers-2 if necessary
182 for(std::size_t i = 1; i <= num_markers - 2; ++i)
184 // offset to desired position
185 float_type d = this->desired_positions[i] - this->actual_positions[i];
187 // offset to next position
188 float_type dp = this->actual_positions[i+1] - this->actual_positions[i];
190 // offset to previous position
191 float_type dm = this->actual_positions[i-1] - this->actual_positions[i];
194 float_type hp = (this->heights[i+1] - this->heights[i]) / dp;
195 float_type hm = (this->heights[i-1] - this->heights[i]) / dm;
197 if((d >= 1 && dp > 1) || (d <= -1 && dm < -1))
199 short sign_d = static_cast<short>(d / std::abs(d));
201 float_type h = this->heights[i] + sign_d / (dp - dm) * ((sign_d - dm)*hp
202 + (dp - sign_d) * hm);
204 // try adjusting heights[i] using p-squared formula
205 if(this->heights[i - 1] < h && h < this->heights[i + 1])
207 this->heights[i] = h;
211 // use linear formula
214 this->heights[i] += hp;
218 this->heights[i] -= hm;
221 this->actual_positions[i] += sign_d;
227 result_type result(dont_care) const
229 // for i in [1,probabilities.size()], return heights[i * 2]
230 detail::times2_iterator idx_begin = detail::make_times2_iterator(1);
231 detail::times2_iterator idx_end = detail::make_times2_iterator(this->probabilities.size() + 1);
234 make_permutation_iterator(this->heights.begin(), idx_begin)
235 , make_permutation_iterator(this->heights.begin(), idx_end)
240 array_type probabilities; // the quantile probabilities
241 array_type heights; // q_i
242 array_type actual_positions; // n_i
243 array_type desired_positions; // d_i
244 array_type positions_increments; // f_i
249 ///////////////////////////////////////////////////////////////////////////////
250 // tag::extended_p_square
254 struct extended_p_square
256 , extended_p_square_probabilities
258 typedef accumulators::impl::extended_p_square_impl<mpl::_1> impl;
260 #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED
261 /// tag::extended_p_square::probabilities named parameter
262 static boost::parameter::keyword<tag::probabilities> const probabilities;
267 ///////////////////////////////////////////////////////////////////////////////
268 // extract::extended_p_square
272 extractor<tag::extended_p_square> const extended_p_square = {};
274 BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square)
277 using extract::extended_p_square;
279 // So that extended_p_square can be automatically substituted with
280 // weighted_extended_p_square when the weight parameter is non-void
282 struct as_weighted_feature<tag::extended_p_square>
284 typedef tag::weighted_extended_p_square type;
288 struct feature_of<tag::weighted_extended_p_square>
289 : feature_of<tag::extended_p_square>
293 }} // namespace boost::accumulators