PageRenderTime 55ms CodeModel.GetById 14ms app.highlight 36ms RepoModel.GetById 2ms app.codeStats 0ms

/Src/Dependencies/Boost/boost/accumulators/statistics/p_square_quantile.hpp

http://hadesmem.googlecode.com/
C++ Header | 257 lines | 171 code | 32 blank | 54 comment | 18 complexity | 3aca4125adf25510506e69b1a5fad761 MD5 | raw file
  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
 26namespace boost { namespace accumulators
 27{
 28
 29namespace 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//
207namespace 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//
229namespace 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    BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_quantile)
235    BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_quantile_for_median)
236}
237
238using extract::p_square_quantile;
239using extract::p_square_quantile_for_median;
240
241// So that p_square_quantile can be automatically substituted with
242// weighted_p_square_quantile when the weight parameter is non-void
243template<>
244struct as_weighted_feature<tag::p_square_quantile>
245{
246    typedef tag::weighted_p_square_quantile type;
247};
248
249template<>
250struct feature_of<tag::weighted_p_square_quantile>
251  : feature_of<tag::p_square_quantile>
252{
253};
254
255}} // namespace boost::accumulators
256
257#endif