boost/accumulators/statistics/p_square_quantile.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 // 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_P_SQUARE_QUANTILE_HPP_DE_01_01_2006
   9 #define BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_QUANTILE_HPP_DE_01_01_2006
  10
  11 #include <cmath>
  12 #include <functional>
  13 #include <boost/array.hpp>
  14 #include <boost/mpl/placeholders.hpp>
  15 #include <boost/type_traits/is_same.hpp>
  16 #include <boost/parameter/keyword.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/framework/depends_on.hpp>
  22 #include <boost/accumulators/statistics_fwd.hpp>
  23 #include <boost/accumulators/statistics/count.hpp>
  24 #include <boost/accumulators/statistics/parameters/quantile_probability.hpp>
  25
  26 namespace boost { namespace accumulators
  27 {
  28
  29 namespace impl
  30 {
  31     ///////////////////////////////////////////////////////////////////////////////
  32     // p_square_quantile_impl
  33     //  single quantile estimation
  34     /**
  35         @brief Single quantile estimation with the \f$P^2\f$ algorithm
  36
  37         The \f$P^2\f$ algorithm estimates a quantile dynamically without storing samples. Instead of
  38         storing the whole sample cumulative distribution, only five points (markers) are stored. The heights
  39         of these markers are the minimum and the maximum of the samples and the current estimates of the
  40         \f$(p/2)\f$-, \f$p\f$- and \f$(1+p)/2\f$-quantiles. Their positions are equal to the number
  41         of samples that are smaller or equal to the markers. Each time a new samples is recorded, the
  42         positions of the markers are updated and if necessary their heights are adjusted using a piecewise-
  43         parabolic formula.
  44
  45         For further details, see
  46
  47         R. Jain and I. Chlamtac, The P^2 algorithmus fordynamic calculation of quantiles and
  48         histograms without storing observations, Communications of the ACM,
  49         Volume 28 (October), Number 10, 1985, p. 1076-1085.
  50
  51         @param quantile_probability
  52     */
  53     template<typename Sample, typename Impl>
  54     struct p_square_quantile_impl
  55       : accumulator_base
  56     {
  57         typedef typename numeric::functional::average<Sample, std::size_t>::result_type float_type;
  58         typedef array<float_type, 5> array_type;
  59         // for boost::result_of
  60         typedef float_type result_type;
  61
  62         template<typename Args>
  63         p_square_quantile_impl(Args const &args)
  64           : p(is_same<Impl, for_median>::value ? 0.5 : args[quantile_probability | 0.5])
  65           , heights()
  66           , actual_positions()
  67           , desired_positions()
  68           , positions_increments()
  69         {
  70             for(std::size_t i = 0; i < 5; ++i)
  71             {
  72                 this->actual_positions[i] = i + 1;
  73             }
  74
  75             this->desired_positions[0] = 1.;
  76             this->desired_positions[1] = 1. + 2. * this->p;
  77             this->desired_positions[2] = 1. + 4. * this->p;
  78             this->desired_positions[3] = 3. + 2. * this->p;
  79             this->desired_positions[4] = 5.;
  80
  81             this->positions_increments[0] = 0.;
  82             this->positions_increments[1] = this->p / 2.;
  83             this->positions_increments[2] = this->p;
  84             this->positions_increments[3] = (1. + this->p) / 2.;
  85             this->positions_increments[4] = 1.;
  86         }
  87
  88         template<typename Args>
  89         void operator ()(Args const &args)
  90         {
  91             std::size_t cnt = count(args);
  92
  93             // accumulate 5 first samples
  94             if(cnt <= 5)
  95             {
  96                 this->heights[cnt - 1] = args[sample];
  97
  98                 // complete the initialization of heights by sorting
  99                 if(cnt == 5)
 100                 {
 101                     std::sort(this->heights.begin(), this->heights.end());
 102                 }
 103             }
 104             else
 105             {
 106                 std::size_t sample_cell = 1; // k
 107
 108                 // find cell k such that heights[k-1] <= args[sample] < heights[k] and ajust extreme values
 109                 if (args[sample] < this->heights[0])
 110                 {
 111                     this->heights[0] = args[sample];
 112                     sample_cell = 1;
 113                 }
 114                 else if (this->heights[4] <= args[sample])
 115                 {
 116                     this->heights[4] = args[sample];
 117                     sample_cell = 4;
 118                 }
 119                 else
 120                 {
 121                     typedef typename array_type::iterator iterator;
 122                     iterator it = std::upper_bound(
 123                         this->heights.begin()
 124                       , this->heights.end()
 125                       , args[sample]
 126                     );
 127
 128                     sample_cell = std::distance(this->heights.begin(), it);
 129                 }
 130
 131                 // update positions of markers above sample_cell
 132                 for(std::size_t i = sample_cell; i < 5; ++i)
 133                 {
 134                     ++this->actual_positions[i];
 135                 }
 136
 137                 // update desired positions of all markers
 138                 for(std::size_t i = 0; i < 5; ++i)
 139                 {
 140                     this->desired_positions[i] += this->positions_increments[i];
 141                 }
 142
 143                 // adjust heights and actual positions of markers 1 to 3 if necessary
 144                 for(std::size_t i = 1; i <= 3; ++i)
 145                 {
 146                     // offset to desired positions
 147                     float_type d = this->desired_positions[i] - this->actual_positions[i];
 148
 149                     // offset to next position
 150                     float_type dp = this->actual_positions[i + 1] - this->actual_positions[i];
 151
 152                     // offset to previous position
 153                     float_type dm = this->actual_positions[i - 1] - this->actual_positions[i];
 154
 155                     // height ds
 156                     float_type hp = (this->heights[i + 1] - this->heights[i]) / dp;
 157                     float_type hm = (this->heights[i - 1] - this->heights[i]) / dm;
 158
 159                     if((d >= 1. && dp > 1.) || (d <= -1. && dm < -1.))
 160                     {
 161                         short sign_d = static_cast<short>(d / std::abs(d));
 162
 163                         // try adjusting heights[i] using p-squared formula
 164                         float_type h = this->heights[i] + sign_d / (dp - dm) * ((sign_d - dm) * hp
 165                                      + (dp - sign_d) * hm);
 166
 167                         if(this->heights[i - 1] < h && h < this->heights[i + 1])
 168                         {
 169                             this->heights[i] = h;
 170                         }
 171                         else
 172                         {
 173                             // use linear formula
 174                             if(d > 0)
 175                             {
 176                                 this->heights[i] += hp;
 177                             }
 178                             if(d < 0)
 179                             {
 180                                 this->heights[i] -= hm;
 181                             }
 182                         }
 183                         this->actual_positions[i] += sign_d;
 184                     }
 185                 }
 186             }
 187         }
 188
 189         result_type result(dont_care) const
 190         {
 191             return this->heights[2];
 192         }
 193
 194     private:
 195         float_type p;                    // the quantile probability p
 196         array_type heights;              // q_i
 197         array_type actual_positions;     // n_i
 198         array_type desired_positions;    // n'_i
 199         array_type positions_increments; // dn'_i
 200     };
 201
 202 } // namespace detail
 203
 204 ///////////////////////////////////////////////////////////////////////////////
 205 // tag::p_square_quantile
 206 //
 207 namespace tag
 208 {
 209     struct p_square_quantile
 210       : depends_on<count>
 211     {
 212         /// INTERNAL ONLY
 213         ///
 214         typedef accumulators::impl::p_square_quantile_impl<mpl::_1, regular> impl;
 215     };
 216     struct p_square_quantile_for_median
 217       : depends_on<count>
 218     {
 219         /// INTERNAL ONLY
 220         ///
 221         typedef accumulators::impl::p_square_quantile_impl<mpl::_1, for_median> impl;
 222     };
 223 }
 224
 225 ///////////////////////////////////////////////////////////////////////////////
 226 // extract::p_square_quantile
 227 // extract::p_square_quantile_for_median
 228 //
 229 namespace extract
 230 {
 231     extractor<tag::p_square_quantile> const p_square_quantile = {};
 232     extractor<tag::p_square_quantile_for_median> const p_square_quantile_for_median = {};
 233 }
 234
 235 using extract::p_square_quantile;
 236 using extract::p_square_quantile_for_median;
 237
 238 // So that p_square_quantile can be automatically substituted with
 239 // weighted_p_square_quantile when the weight parameter is non-void
 240 template<>
 241 struct as_weighted_feature<tag::p_square_quantile>
 242 {
 243     typedef tag::weighted_p_square_quantile type;
 244 };
 245
 246 template<>
 247 struct feature_of<tag::weighted_p_square_quantile>
 248   : feature_of<tag::p_square_quantile>
 249 {
 250 };
 251
 252 }} // namespace boost::accumulators
 253
 254 #endif