Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // p_square_quantile.hpp
Chris@16: //
Chris@16: //  Copyright 2005 Daniel Egloff. Distributed under the Boost
Chris@16: //  Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_QUANTILE_HPP_DE_01_01_2006
Chris@16: #define BOOST_ACCUMULATORS_STATISTICS_P_SQUARE_QUANTILE_HPP_DE_01_01_2006
Chris@16: 
Chris@16: #include <cmath>
Chris@16: #include <functional>
Chris@16: #include <boost/array.hpp>
Chris@16: #include <boost/mpl/placeholders.hpp>
Chris@16: #include <boost/type_traits/is_same.hpp>
Chris@16: #include <boost/parameter/keyword.hpp>
Chris@16: #include <boost/accumulators/framework/accumulator_base.hpp>
Chris@16: #include <boost/accumulators/framework/extractor.hpp>
Chris@16: #include <boost/accumulators/numeric/functional.hpp>
Chris@16: #include <boost/accumulators/framework/parameters/sample.hpp>
Chris@16: #include <boost/accumulators/framework/depends_on.hpp>
Chris@16: #include <boost/accumulators/statistics_fwd.hpp>
Chris@16: #include <boost/accumulators/statistics/count.hpp>
Chris@16: #include <boost/accumulators/statistics/parameters/quantile_probability.hpp>
Chris@16: 
Chris@16: namespace boost { namespace accumulators
Chris@16: {
Chris@16: 
Chris@16: namespace impl
Chris@16: {
Chris@16:     ///////////////////////////////////////////////////////////////////////////////
Chris@16:     // p_square_quantile_impl
Chris@16:     //  single quantile estimation
Chris@16:     /**
Chris@16:         @brief Single quantile estimation with the \f$P^2\f$ algorithm
Chris@16: 
Chris@16:         The \f$P^2\f$ algorithm estimates a quantile dynamically without storing samples. Instead of
Chris@16:         storing the whole sample cumulative distribution, only five points (markers) are stored. The heights
Chris@16:         of these markers are the minimum and the maximum of the samples and the current estimates of the
Chris@16:         \f$(p/2)\f$-, \f$p\f$- and \f$(1+p)/2\f$-quantiles. Their positions are equal to the number
Chris@16:         of samples that are smaller or equal to the markers. Each time a new samples is recorded, the
Chris@16:         positions of the markers are updated and if necessary their heights are adjusted using a piecewise-
Chris@16:         parabolic formula.
Chris@16: 
Chris@16:         For further details, see
Chris@16: 
Chris@16:         R. Jain and I. Chlamtac, The P^2 algorithm for dynamic calculation of quantiles and
Chris@16:         histograms without storing observations, Communications of the ACM,
Chris@16:         Volume 28 (October), Number 10, 1985, p. 1076-1085.
Chris@16: 
Chris@16:         @param quantile_probability
Chris@16:     */
Chris@16:     template<typename Sample, typename Impl>
Chris@16:     struct p_square_quantile_impl
Chris@16:       : accumulator_base
Chris@16:     {
Chris@16:         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
Chris@16:         typedef array<float_type, 5> array_type;
Chris@16:         // for boost::result_of
Chris@16:         typedef float_type result_type;
Chris@16: 
Chris@16:         template<typename Args>
Chris@16:         p_square_quantile_impl(Args const &args)
Chris@16:           : p(is_same<Impl, for_median>::value ? 0.5 : args[quantile_probability | 0.5])
Chris@16:           , heights()
Chris@16:           , actual_positions()
Chris@16:           , desired_positions()
Chris@16:           , positions_increments()
Chris@16:         {
Chris@16:             for(std::size_t i = 0; i < 5; ++i)
Chris@16:             {
Chris@16:                 this->actual_positions[i] = i + 1.;
Chris@16:             }
Chris@16: 
Chris@16:             this->desired_positions[0] = 1.;
Chris@16:             this->desired_positions[1] = 1. + 2. * this->p;
Chris@16:             this->desired_positions[2] = 1. + 4. * this->p;
Chris@16:             this->desired_positions[3] = 3. + 2. * this->p;
Chris@16:             this->desired_positions[4] = 5.;
Chris@16: 
Chris@16:             this->positions_increments[0] = 0.;
Chris@16:             this->positions_increments[1] = this->p / 2.;
Chris@16:             this->positions_increments[2] = this->p;
Chris@16:             this->positions_increments[3] = (1. + this->p) / 2.;
Chris@16:             this->positions_increments[4] = 1.;
Chris@16:         }
Chris@16: 
Chris@16:         template<typename Args>
Chris@16:         void operator ()(Args const &args)
Chris@16:         {
Chris@16:             std::size_t cnt = count(args);
Chris@16: 
Chris@16:             // accumulate 5 first samples
Chris@16:             if(cnt <= 5)
Chris@16:             {
Chris@16:                 this->heights[cnt - 1] = args[sample];
Chris@16: 
Chris@16:                 // complete the initialization of heights by sorting
Chris@16:                 if(cnt == 5)
Chris@16:                 {
Chris@16:                     std::sort(this->heights.begin(), this->heights.end());
Chris@16:                 }
Chris@16:             }
Chris@16:             else
Chris@16:             {
Chris@16:                 std::size_t sample_cell = 1; // k
Chris@16: 
Chris@16:                 // find cell k such that heights[k-1] <= args[sample] < heights[k] and adjust extreme values
Chris@16:                 if (args[sample] < this->heights[0])
Chris@16:                 {
Chris@16:                     this->heights[0] = args[sample];
Chris@16:                     sample_cell = 1;
Chris@16:                 }
Chris@16:                 else if (this->heights[4] <= args[sample])
Chris@16:                 {
Chris@16:                     this->heights[4] = args[sample];
Chris@16:                     sample_cell = 4;
Chris@16:                 }
Chris@16:                 else
Chris@16:                 {
Chris@16:                     typedef typename array_type::iterator iterator;
Chris@16:                     iterator it = std::upper_bound(
Chris@16:                         this->heights.begin()
Chris@16:                       , this->heights.end()
Chris@16:                       , args[sample]
Chris@16:                     );
Chris@16: 
Chris@16:                     sample_cell = std::distance(this->heights.begin(), it);
Chris@16:                 }
Chris@16: 
Chris@16:                 // update positions of markers above sample_cell
Chris@16:                 for(std::size_t i = sample_cell; i < 5; ++i)
Chris@16:                 {
Chris@16:                     ++this->actual_positions[i];
Chris@16:                 }
Chris@16: 
Chris@16:                 // update desired positions of all markers
Chris@16:                 for(std::size_t i = 0; i < 5; ++i)
Chris@16:                 {
Chris@16:                     this->desired_positions[i] += this->positions_increments[i];
Chris@16:                 }
Chris@16: 
Chris@16:                 // adjust heights and actual positions of markers 1 to 3 if necessary
Chris@16:                 for(std::size_t i = 1; i <= 3; ++i)
Chris@16:                 {
Chris@16:                     // offset to desired positions
Chris@16:                     float_type d = this->desired_positions[i] - this->actual_positions[i];
Chris@16: 
Chris@16:                     // offset to next position
Chris@16:                     float_type dp = this->actual_positions[i + 1] - this->actual_positions[i];
Chris@16: 
Chris@16:                     // offset to previous position
Chris@16:                     float_type dm = this->actual_positions[i - 1] - this->actual_positions[i];
Chris@16: 
Chris@16:                     // height ds
Chris@16:                     float_type hp = (this->heights[i + 1] - this->heights[i]) / dp;
Chris@16:                     float_type hm = (this->heights[i - 1] - this->heights[i]) / dm;
Chris@16: 
Chris@16:                     if((d >= 1. && dp > 1.) || (d <= -1. && dm < -1.))
Chris@16:                     {
Chris@16:                         short sign_d = static_cast<short>(d / std::abs(d));
Chris@16: 
Chris@16:                         // try adjusting heights[i] using p-squared formula
Chris@16:                         float_type h = this->heights[i] + sign_d / (dp - dm) * ((sign_d - dm) * hp
Chris@16:                                      + (dp - sign_d) * hm);
Chris@16: 
Chris@16:                         if(this->heights[i - 1] < h && h < this->heights[i + 1])
Chris@16:                         {
Chris@16:                             this->heights[i] = h;
Chris@16:                         }
Chris@16:                         else
Chris@16:                         {
Chris@16:                             // use linear formula
Chris@16:                             if(d > 0)
Chris@16:                             {
Chris@16:                                 this->heights[i] += hp;
Chris@16:                             }
Chris@16:                             if(d < 0)
Chris@16:                             {
Chris@16:                                 this->heights[i] -= hm;
Chris@16:                             }
Chris@16:                         }
Chris@16:                         this->actual_positions[i] += sign_d;
Chris@16:                     }
Chris@16:                 }
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         result_type result(dont_care) const
Chris@16:         {
Chris@16:             return this->heights[2];
Chris@16:         }
Chris@16: 
Chris@16:     private:
Chris@16:         float_type p;                    // the quantile probability p
Chris@16:         array_type heights;              // q_i
Chris@16:         array_type actual_positions;     // n_i
Chris@16:         array_type desired_positions;    // n'_i
Chris@16:         array_type positions_increments; // dn'_i
Chris@16:     };
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // tag::p_square_quantile
Chris@16: //
Chris@16: namespace tag
Chris@16: {
Chris@16:     struct p_square_quantile
Chris@16:       : depends_on<count>
Chris@16:     {
Chris@16:         /// INTERNAL ONLY
Chris@16:         ///
Chris@16:         typedef accumulators::impl::p_square_quantile_impl<mpl::_1, regular> impl;
Chris@16:     };
Chris@16:     struct p_square_quantile_for_median
Chris@16:       : depends_on<count>
Chris@16:     {
Chris@16:         /// INTERNAL ONLY
Chris@16:         ///
Chris@16:         typedef accumulators::impl::p_square_quantile_impl<mpl::_1, for_median> impl;
Chris@16:     };
Chris@16: }
Chris@16: 
Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // extract::p_square_quantile
Chris@16: // extract::p_square_quantile_for_median
Chris@16: //
Chris@16: namespace extract
Chris@16: {
Chris@16:     extractor<tag::p_square_quantile> const p_square_quantile = {};
Chris@16:     extractor<tag::p_square_quantile_for_median> const p_square_quantile_for_median = {};
Chris@16: 
Chris@16:     BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_quantile)
Chris@16:     BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_quantile_for_median)
Chris@16: }
Chris@16: 
Chris@16: using extract::p_square_quantile;
Chris@16: using extract::p_square_quantile_for_median;
Chris@16: 
Chris@16: // So that p_square_quantile can be automatically substituted with
Chris@16: // weighted_p_square_quantile when the weight parameter is non-void
Chris@16: template<>
Chris@16: struct as_weighted_feature<tag::p_square_quantile>
Chris@16: {
Chris@16:     typedef tag::weighted_p_square_quantile type;
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct feature_of<tag::weighted_p_square_quantile>
Chris@16:   : feature_of<tag::p_square_quantile>
Chris@16: {
Chris@16: };
Chris@16: 
Chris@16: }} // namespace boost::accumulators
Chris@16: 
Chris@16: #endif