PageRenderTime 45ms CodeModel.GetById 20ms app.highlight 21ms RepoModel.GetById 1ms app.codeStats 0ms

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

http://hadesmem.googlecode.com/
C++ Header | 293 lines | 195 code | 38 blank | 60 comment | 21 complexity | cd3bb37c768c007548bdfc8e8e30e6a6 MD5 | raw file
  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
 30namespace boost { namespace accumulators
 31{
 32///////////////////////////////////////////////////////////////////////////////
 33// probabilities named parameter
 34//
 35BOOST_PARAMETER_NESTED_KEYWORD(tag, extended_p_square_probabilities, probabilities)
 36
 37namespace 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//
250namespace 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//
268namespace extract
269{
270    extractor<tag::extended_p_square> const extended_p_square = {};
271
272    BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square)
273}
274
275using extract::extended_p_square;
276
277// So that extended_p_square can be automatically substituted with
278// weighted_extended_p_square when the weight parameter is non-void
279template<>
280struct as_weighted_feature<tag::extended_p_square>
281{
282    typedef tag::weighted_extended_p_square type;
283};
284
285template<>
286struct feature_of<tag::weighted_extended_p_square>
287  : feature_of<tag::extended_p_square>
288{
289};
290
291}} // namespace boost::accumulators
292
293#endif