Statistics
| Revision:

root / tmp / org.txm.statsengine.r.core.win32 / res / win32 / library / BH / include / boost / accumulators / statistics / p_square_quantile.hpp @ 2486

History | View | Annotate | Download (9.1 kB)

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 algorithm for dynamic 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::fdiv<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 adjust 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
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_quantile)
235
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(p_square_quantile_for_median)
236
}
237

    
238
using extract::p_square_quantile;
239
using 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
243
template<>
244
struct as_weighted_feature<tag::p_square_quantile>
245
{
246
    typedef tag::weighted_p_square_quantile type;
247
};
248

    
249
template<>
250
struct feature_of<tag::weighted_p_square_quantile>
251
  : feature_of<tag::p_square_quantile>
252
{
253
};
254

    
255
}} // namespace boost::accumulators
256

    
257
#endif