1 ///////////////////////////////////////////////////////////////////////////////
2 // p_square_cumulative_distribution.hpp
4 // Copyright 2005 Daniel Egloff, Olivier Gygi. 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_P_SQUARE_CUMUL_DIST_HPP_DE_01_01_2006
9 #define BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_CUMUL_DIST_HPP_DE_01_01_2006
13 #include <boost/parameter/keyword.hpp>
14 #include <boost/range.hpp>
15 #include <boost/mpl/placeholders.hpp>
16 #include <boost/accumulators/accumulators_fwd.hpp>
17 #include <boost/accumulators/framework/accumulator_base.hpp>
18 #include <boost/accumulators/framework/extractor.hpp>
19 #include <boost/accumulators/numeric/functional.hpp>
20 #include <boost/accumulators/framework/parameters/sample.hpp>
21 #include <boost/accumulators/statistics_fwd.hpp>
22 #include <boost/accumulators/statistics/count.hpp>
24 namespace boost { namespace accumulators
26 ///////////////////////////////////////////////////////////////////////////////
27 // num_cells named parameter
29 BOOST_PARAMETER_NESTED_KEYWORD(tag, p_square_cumulative_distribution_num_cells, num_cells)
31 BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_cumulative_distribution_num_cells)
35 ///////////////////////////////////////////////////////////////////////////////
36 // p_square_cumulative_distribution_impl
37 // cumulative_distribution calculation (as histogram)
39 @brief Histogram calculation of the cumulative distribution with the \f$P^2\f$ algorithm
41 A histogram of the sample cumulative distribution is computed dynamically without storing samples
42 based on the \f$ P^2 \f$ algorithm. The returned histogram has a specifiable amount (num_cells)
43 equiprobable (and not equal-sized) cells.
45 For further details, see
47 R. Jain and I. Chlamtac, The P^2 algorithm for dynamic calculation of quantiles and
48 histograms without storing observations, Communications of the ACM,
49 Volume 28 (October), Number 10, 1985, p. 1076-1085.
51 @param p_square_cumulative_distribution_num_cells.
53 template<typename Sample>
54 struct p_square_cumulative_distribution_impl
57 typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
58 typedef std::vector<float_type> array_type;
59 typedef std::vector<std::pair<float_type, float_type> > histogram_type;
60 // for boost::result_of
61 typedef iterator_range<typename histogram_type::iterator> result_type;
63 template<typename Args>
64 p_square_cumulative_distribution_impl(Args const &args)
65 : num_cells(args[p_square_cumulative_distribution_num_cells])
66 , heights(num_cells + 1)
67 , actual_positions(num_cells + 1)
68 , desired_positions(num_cells + 1)
69 , positions_increments(num_cells + 1)
70 , histogram(num_cells + 1)
73 std::size_t b = this->num_cells;
75 for (std::size_t i = 0; i < b + 1; ++i)
77 this->actual_positions[i] = i + 1.;
78 this->desired_positions[i] = i + 1.;
79 this->positions_increments[i] = numeric::fdiv(i, b);
83 template<typename Args>
84 void operator ()(Args const &args)
86 this->is_dirty = true;
88 std::size_t cnt = count(args);
89 std::size_t sample_cell = 1; // k
90 std::size_t b = this->num_cells;
92 // accumulate num_cells + 1 first samples
95 this->heights[cnt - 1] = args[sample];
97 // complete the initialization of heights by sorting
100 std::sort(this->heights.begin(), this->heights.end());
105 // find cell k such that heights[k-1] <= args[sample] < heights[k] and adjust extreme values
106 if (args[sample] < this->heights[0])
108 this->heights[0] = args[sample];
111 else if (this->heights[b] <= args[sample])
113 this->heights[b] = args[sample];
118 typename array_type::iterator it;
119 it = std::upper_bound(
120 this->heights.begin()
121 , this->heights.end()
125 sample_cell = std::distance(this->heights.begin(), it);
128 // increment positions of markers above sample_cell
129 for (std::size_t i = sample_cell; i < b + 1; ++i)
131 ++this->actual_positions[i];
134 // update desired position of markers 2 to num_cells + 1
135 // (desired position of first marker is always 1)
136 for (std::size_t i = 1; i < b + 1; ++i)
138 this->desired_positions[i] += this->positions_increments[i];
141 // adjust heights of markers 2 to num_cells if necessary
142 for (std::size_t i = 1; i < b; ++i)
144 // offset to desire position
145 float_type d = this->desired_positions[i] - this->actual_positions[i];
147 // offset to next position
148 float_type dp = this->actual_positions[i + 1] - this->actual_positions[i];
150 // offset to previous position
151 float_type dm = this->actual_positions[i - 1] - this->actual_positions[i];
154 float_type hp = (this->heights[i + 1] - this->heights[i]) / dp;
155 float_type hm = (this->heights[i - 1] - this->heights[i]) / dm;
157 if ( ( d >= 1. && dp > 1. ) || ( d <= -1. && dm < -1. ) )
159 short sign_d = static_cast<short>(d / std::abs(d));
161 // try adjusting heights[i] using p-squared formula
162 float_type h = this->heights[i] + sign_d / (dp - dm) * ( (sign_d - dm) * hp + (dp - sign_d) * hm );
164 if ( this->heights[i - 1] < h && h < this->heights[i + 1] )
166 this->heights[i] = h;
170 // use linear formula
173 this->heights[i] += hp;
177 this->heights[i] -= hm;
180 this->actual_positions[i] += sign_d;
186 template<typename Args>
187 result_type result(Args const &args) const
191 this->is_dirty = false;
193 // creates a vector of std::pair where each pair i holds
194 // the values heights[i] (x-axis of histogram) and
195 // actual_positions[i] / cnt (y-axis of histogram)
197 std::size_t cnt = count(args);
199 for (std::size_t i = 0; i < this->histogram.size(); ++i)
201 this->histogram[i] = std::make_pair(this->heights[i], numeric::fdiv(this->actual_positions[i], cnt));
205 return make_iterator_range(this->histogram);
209 std::size_t num_cells; // number of cells b
210 array_type heights; // q_i
211 array_type actual_positions; // n_i
212 array_type desired_positions; // n'_i
213 array_type positions_increments; // dn'_i
214 mutable histogram_type histogram; // histogram
215 mutable bool is_dirty;
218 } // namespace detail
220 ///////////////////////////////////////////////////////////////////////////////
221 // tag::p_square_cumulative_distribution
225 struct p_square_cumulative_distribution
227 , p_square_cumulative_distribution_num_cells
231 typedef accumulators::impl::p_square_cumulative_distribution_impl<mpl::_1> impl;
235 ///////////////////////////////////////////////////////////////////////////////
236 // extract::p_square_cumulative_distribution
240 extractor<tag::p_square_cumulative_distribution> const p_square_cumulative_distribution = {};
242 BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_cumulative_distribution)
245 using extract::p_square_cumulative_distribution;
247 // So that p_square_cumulative_distribution can be automatically substituted with
248 // weighted_p_square_cumulative_distribution when the weight parameter is non-void
250 struct as_weighted_feature<tag::p_square_cumulative_distribution>
252 typedef tag::weighted_p_square_cumulative_distribution type;
256 struct feature_of<tag::weighted_p_square_cumulative_distribution>
257 : feature_of<tag::p_square_cumulative_distribution>
261 }} // namespace boost::accumulators