boost/accumulators/statistics/extended_p_square.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 // extended_p_square.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_HPP_DE_01_01_2006
   9 #define BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_HPP_DE_01_01_2006
  10
  11 #include <vector>
  12 #include <functional>
  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/framework/accumulator_base.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>
  29
  30 namespace boost { namespace accumulators
  31 {
  32 ///////////////////////////////////////////////////////////////////////////////
  33 // probabilities named parameter
  34 //
  35 BOOST_PARAMETER_NESTED_KEYWORD(tag, extended_p_square_probabilities, probabilities)
  36
  37 namespace impl
  38 {
  39     ///////////////////////////////////////////////////////////////////////////////
  40     // extended_p_square_impl
  41     //  multiple quantile estimation
  42     /**
  43         @brief Multiple quantile estimation with the extended \f$P^2\f$ algorithm
  44
  45         Extended \f$P^2\f$ algorithm for estimation of several quantiles without storing samples.
  46         Assume that \f$m\f$ quantiles \f$\xi_{p_1}, \ldots, \xi_{p_m}\f$ are to be estimated.
  47         Instead of storing the whole sample cumulative distribution, the algorithm maintains only
  48         \f$m+2\f$ principal markers and \f$m+1\f$ middle markers, whose positions are updated
  49         with each sample and whose heights are adjusted (if necessary) using a piecewise-parablic
  50         formula. The heights of these central markers are the current estimates of the quantiles
  51         and returned as an iterator range.
  52
  53         For further details, see
  54
  55         K. E. E. Raatikainen, Simultaneous estimation of several quantiles, Simulation, Volume 49,
  56         Number 4 (October), 1986, p. 159-164.
  57
  58         The extended \f$ P^2 \f$ algorithm generalizess the \f$ P^2 \f$ algorithm of
  59
  60         R. Jain and I. Chlamtac, The P^2 algorithmus for dynamic calculation of quantiles and
  61         histograms without storing observations, Communications of the ACM,
  62         Volume 28 (October), Number 10, 1985, p. 1076-1085.
  63
  64         @param extended_p_square_probabilities A vector of quantile probabilities.
  65     */
  66     template<typename Sample>
  67     struct extended_p_square_impl
  68       : accumulator_base
  69     {
  70         typedef typename numeric::functional::average<Sample, std::size_t>::result_type float_type;
  71         typedef std::vector<float_type> array_type;
  72         // for boost::result_of
  73         typedef iterator_range<
  74             detail::lvalue_index_iterator<
  75                 permutation_iterator<
  76                     typename array_type::const_iterator
  77                   , detail::times2_iterator
  78                 >
  79             >
  80         > result_type;
  81
  82         template<typename Args>
  83         extended_p_square_impl(Args const &args)
  84           : probabilities(
  85                 boost::begin(args[extended_p_square_probabilities])
  86               , boost::end(args[extended_p_square_probabilities])
  87             )
  88           , heights(2 * probabilities.size() + 3)
  89           , actual_positions(heights.size())
  90           , desired_positions(heights.size())
  91           , positions_increments(heights.size())
  92         {
  93             std::size_t num_quantiles = this->probabilities.size();
  94             std::size_t num_markers = this->heights.size();
  95
  96             for(std::size_t i = 0; i < num_markers; ++i)
  97             {
  98                 this->actual_positions[i] = i + 1;
  99             }
 100
 101             this->positions_increments[0] = 0.;
 102             this->positions_increments[num_markers - 1] = 1.;
 103
 104             for(std::size_t i = 0; i < num_quantiles; ++i)
 105             {
 106                 this->positions_increments[2 * i + 2] = probabilities[i];
 107             }
 108
 109             for(std::size_t i = 0; i <= num_quantiles; ++i)
 110             {
 111                 this->positions_increments[2 * i + 1] =
 112                     0.5 * (this->positions_increments[2 * i] + this->positions_increments[2 * i + 2]);
 113             }
 114
 115             for(std::size_t i = 0; i < num_markers; ++i)
 116             {
 117                 this->desired_positions[i] = 1. + 2. * (num_quantiles + 1.) * this->positions_increments[i];
 118             }
 119         }
 120
 121         template<typename Args>
 122         void operator ()(Args const &args)
 123         {
 124             std::size_t cnt = count(args);
 125
 126             // m+2 principal markers and m+1 middle markers
 127             std::size_t num_markers = 2 * this->probabilities.size() + 3;
 128
 129             // first accumulate num_markers samples
 130             if(cnt <= num_markers)
 131             {
 132                 this->heights[cnt - 1] = args[sample];
 133
 134                 // complete the initialization of heights by sorting
 135                 if(cnt == num_markers)
 136                 {
 137                     std::sort(this->heights.begin(), this->heights.end());
 138                 }
 139             }
 140             else
 141             {
 142                 std::size_t sample_cell = 1;
 143
 144                 // find cell k = sample_cell such that heights[k-1] <= sample < heights[k]
 145                 if(args[sample] < this->heights[0])
 146                 {
 147                     this->heights[0] = args[sample];
 148                     sample_cell = 1;
 149                 }
 150                 else if(args[sample] >= this->heights[num_markers - 1])
 151                 {
 152                     this->heights[num_markers - 1] = args[sample];
 153                     sample_cell = num_markers - 1;
 154                 }
 155                 else
 156                 {
 157                     typedef typename array_type::iterator iterator;
 158                     iterator it = std::upper_bound(
 159                         this->heights.begin()
 160                       , this->heights.end()
 161                       , args[sample]
 162                     );
 163
 164                     sample_cell = std::distance(this->heights.begin(), it);
 165                 }
 166
 167                 // update actual positions of all markers above sample_cell index
 168                 for(std::size_t i = sample_cell; i < num_markers; ++i)
 169                 {
 170                     ++this->actual_positions[i];
 171                 }
 172
 173                 // update desired positions of all markers
 174                 for(std::size_t i = 0; i < num_markers; ++i)
 175                 {
 176                     this->desired_positions[i] += this->positions_increments[i];
 177                 }
 178
 179                 // adjust heights and actual positions of markers 1 to num_markers-2 if necessary
 180                 for(std::size_t i = 1; i <= num_markers - 2; ++i)
 181                 {
 182                     // offset to desired position
 183                     float_type d = this->desired_positions[i] - this->actual_positions[i];
 184
 185                     // offset to next position
 186                     float_type dp = this->actual_positions[i+1] - this->actual_positions[i];
 187
 188                     // offset to previous position
 189                     float_type dm = this->actual_positions[i-1] - this->actual_positions[i];
 190
 191                     // height ds
 192                     float_type hp = (this->heights[i+1] - this->heights[i]) / dp;
 193                     float_type hm = (this->heights[i-1] - this->heights[i]) / dm;
 194
 195                     if((d >= 1 && dp > 1) || (d <= -1 && dm < -1))
 196                     {
 197                         short sign_d = static_cast<short>(d / std::abs(d));
 198
 199                         float_type h = this->heights[i] + sign_d / (dp - dm) * ((sign_d - dm)*hp
 200                                      + (dp - sign_d) * hm);
 201
 202                         // try adjusting heights[i] using p-squared formula
 203                         if(this->heights[i - 1] < h && h < this->heights[i + 1])
 204                         {
 205                             this->heights[i] = h;
 206                         }
 207                         else
 208                         {
 209                             // use linear formula
 210                             if(d > 0)
 211                             {
 212                                 this->heights[i] += hp;
 213                             }
 214                             if(d < 0)
 215                             {
 216                                 this->heights[i] -= hm;
 217                             }
 218                         }
 219                         this->actual_positions[i] += sign_d;
 220                     }
 221                 }
 222             }
 223         }
 224
 225         result_type result(dont_care) const
 226         {
 227             // for i in [1,probabilities.size()], return heights[i * 2]
 228             detail::times2_iterator idx_begin = detail::make_times2_iterator(1);
 229             detail::times2_iterator idx_end = detail::make_times2_iterator(this->probabilities.size() + 1);
 230
 231             return result_type(
 232                 make_permutation_iterator(this->heights.begin(), idx_begin)
 233               , make_permutation_iterator(this->heights.begin(), idx_end)
 234             );
 235         }
 236
 237     private:
 238         array_type probabilities;         // the quantile probabilities
 239         array_type heights;               // q_i
 240         array_type actual_positions;      // n_i
 241         array_type desired_positions;     // d_i
 242         array_type positions_increments;  // f_i
 243     };
 244
 245 } // namespace impl
 246
 247 ///////////////////////////////////////////////////////////////////////////////
 248 // tag::extended_p_square
 249 //
 250 namespace tag
 251 {
 252     struct extended_p_square
 253       : depends_on<count>
 254       , extended_p_square_probabilities
 255     {
 256         typedef accumulators::impl::extended_p_square_impl<mpl::_1> impl;
 257
 258         #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED
 259         /// tag::extended_p_square::probabilities named paramter
 260         static boost::parameter::keyword<tag::probabilities> const probabilities;
 261         #endif
 262     };
 263 }
 264
 265 ///////////////////////////////////////////////////////////////////////////////
 266 // extract::extended_p_square
 267 //
 268 namespace extract
 269 {
 270     extractor<tag::extended_p_square> const extended_p_square = {};
 271 }
 272
 273 using extract::extended_p_square;
 274
 275 // So that extended_p_square can be automatically substituted with
 276 // weighted_extended_p_square when the weight parameter is non-void
 277 template<>
 278 struct as_weighted_feature<tag::extended_p_square>
 279 {
 280     typedef tag::weighted_extended_p_square type;
 281 };
 282
 283 template<>
 284 struct feature_of<tag::weighted_extended_p_square>
 285   : feature_of<tag::extended_p_square>
 286 {
 287 };
 288
 289 }} // namespace boost::accumulators
 290
 291 #endif