PageRenderTime 35ms CodeModel.GetById 15ms app.highlight 16ms RepoModel.GetById 1ms app.codeStats 0ms

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

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