Statistics
| Revision:

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

History | View | Annotate | Download (10.5 kB)

1 2486 sjacqu01
///////////////////////////////////////////////////////////////////////////////
2 2486 sjacqu01
// weighted_p_square_quantile.hpp
3 2486 sjacqu01
//
4 2486 sjacqu01
//  Copyright 2005 Daniel Egloff. Distributed under the Boost
5 2486 sjacqu01
//  Software License, Version 1.0. (See accompanying file
6 2486 sjacqu01
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
7 2486 sjacqu01
8 2486 sjacqu01
#ifndef BOOST_ACCUMULATORS_STATISTICS_WEIGHTED_P_SQUARE_QUANTILE_HPP_DE_01_01_2006
9 2486 sjacqu01
#define BOOST_ACCUMULATORS_STATISTICS_WEIGHTED_P_SQUARE_QUANTILE_HPP_DE_01_01_2006
10 2486 sjacqu01
11 2486 sjacqu01
#include <cmath>
12 2486 sjacqu01
#include <functional>
13 2486 sjacqu01
#include <boost/array.hpp>
14 2486 sjacqu01
#include <boost/parameter/keyword.hpp>
15 2486 sjacqu01
#include <boost/mpl/placeholders.hpp>
16 2486 sjacqu01
#include <boost/type_traits/is_same.hpp>
17 2486 sjacqu01
#include <boost/accumulators/framework/accumulator_base.hpp>
18 2486 sjacqu01
#include <boost/accumulators/framework/extractor.hpp>
19 2486 sjacqu01
#include <boost/accumulators/numeric/functional.hpp>
20 2486 sjacqu01
#include <boost/accumulators/framework/parameters/sample.hpp>
21 2486 sjacqu01
#include <boost/accumulators/statistics_fwd.hpp>
22 2486 sjacqu01
#include <boost/accumulators/statistics/count.hpp>
23 2486 sjacqu01
#include <boost/accumulators/statistics/sum.hpp>
24 2486 sjacqu01
#include <boost/accumulators/statistics/parameters/quantile_probability.hpp>
25 2486 sjacqu01
26 2486 sjacqu01
namespace boost { namespace accumulators
27 2486 sjacqu01
{
28 2486 sjacqu01
29 2486 sjacqu01
namespace impl {
30 2486 sjacqu01
    ///////////////////////////////////////////////////////////////////////////////
31 2486 sjacqu01
    // weighted_p_square_quantile_impl
32 2486 sjacqu01
    //  single quantile estimation with weighted samples
33 2486 sjacqu01
    /**
34 2486 sjacqu01
        @brief Single quantile estimation with the \f$P^2\f$ algorithm for weighted samples
35 2486 sjacqu01

36 2486 sjacqu01
        This version of the \f$P^2\f$ algorithm extends the \f$P^2\f$ algorithm to support weighted samples.
37 2486 sjacqu01
        The \f$P^2\f$ algorithm estimates a quantile dynamically without storing samples. Instead of
38 2486 sjacqu01
        storing the whole sample cumulative distribution, only five points (markers) are stored. The heights
39 2486 sjacqu01
        of these markers are the minimum and the maximum of the samples and the current estimates of the
40 2486 sjacqu01
        \f$(p/2)\f$-, \f$p\f$ - and \f$(1+p)/2\f$ -quantiles. Their positions are equal to the number
41 2486 sjacqu01
        of samples that are smaller or equal to the markers. Each time a new sample is added, the
42 2486 sjacqu01
        positions of the markers are updated and if necessary their heights are adjusted using a piecewise-
43 2486 sjacqu01
        parabolic formula.
44 2486 sjacqu01

45 2486 sjacqu01
        For further details, see
46 2486 sjacqu01

47 2486 sjacqu01
        R. Jain and I. Chlamtac, The P^2 algorithm for dynamic calculation of quantiles and
48 2486 sjacqu01
        histograms without storing observations, Communications of the ACM,
49 2486 sjacqu01
        Volume 28 (October), Number 10, 1985, p. 1076-1085.
50 2486 sjacqu01

51 2486 sjacqu01
        @param quantile_probability
52 2486 sjacqu01
    */
53 2486 sjacqu01
    template<typename Sample, typename Weight, typename Impl>
54 2486 sjacqu01
    struct weighted_p_square_quantile_impl
55 2486 sjacqu01
      : accumulator_base
56 2486 sjacqu01
    {
57 2486 sjacqu01
        typedef typename numeric::functional::multiplies<Sample, Weight>::result_type weighted_sample;
58 2486 sjacqu01
        typedef typename numeric::functional::fdiv<weighted_sample, std::size_t>::result_type float_type;
59 2486 sjacqu01
        typedef array<float_type, 5> array_type;
60 2486 sjacqu01
        // for boost::result_of
61 2486 sjacqu01
        typedef float_type result_type;
62 2486 sjacqu01
63 2486 sjacqu01
        template<typename Args>
64 2486 sjacqu01
        weighted_p_square_quantile_impl(Args const &args)
65 2486 sjacqu01
          : p(is_same<Impl, for_median>::value ? 0.5 : args[quantile_probability | 0.5])
66 2486 sjacqu01
          , heights()
67 2486 sjacqu01
          , actual_positions()
68 2486 sjacqu01
          , desired_positions()
69 2486 sjacqu01
        {
70 2486 sjacqu01
        }
71 2486 sjacqu01
72 2486 sjacqu01
        template<typename Args>
73 2486 sjacqu01
        void operator ()(Args const &args)
74 2486 sjacqu01
        {
75 2486 sjacqu01
            std::size_t cnt = count(args);
76 2486 sjacqu01
77 2486 sjacqu01
            // accumulate 5 first samples
78 2486 sjacqu01
            if (cnt <= 5)
79 2486 sjacqu01
            {
80 2486 sjacqu01
                this->heights[cnt - 1] = args[sample];
81 2486 sjacqu01
82 2486 sjacqu01
                // In this initialization phase, actual_positions stores the weights of the
83 2486 sjacqu01
                // initial samples that are needed at the end of the initialization phase to
84 2486 sjacqu01
                // compute the correct initial positions of the markers.
85 2486 sjacqu01
                this->actual_positions[cnt - 1] = args[weight];
86 2486 sjacqu01
87 2486 sjacqu01
                // complete the initialization of heights and actual_positions by sorting
88 2486 sjacqu01
                if (cnt == 5)
89 2486 sjacqu01
                {
90 2486 sjacqu01
                    // TODO: we need to sort the initial samples (in heights) in ascending order and
91 2486 sjacqu01
                    // sort their weights (in actual_positions) the same way. The following lines do
92 2486 sjacqu01
                    // it, but there must be a better and more efficient way of doing this.
93 2486 sjacqu01
                    typename array_type::iterator it_begin, it_end, it_min;
94 2486 sjacqu01
95 2486 sjacqu01
                    it_begin = this->heights.begin();
96 2486 sjacqu01
                    it_end   = this->heights.end();
97 2486 sjacqu01
98 2486 sjacqu01
                    std::size_t pos = 0;
99 2486 sjacqu01
100 2486 sjacqu01
                    while (it_begin != it_end)
101 2486 sjacqu01
                    {
102 2486 sjacqu01
                        it_min = std::min_element(it_begin, it_end);
103 2486 sjacqu01
                        std::size_t d = std::distance(it_begin, it_min);
104 2486 sjacqu01
                        std::swap(*it_begin, *it_min);
105 2486 sjacqu01
                        std::swap(this->actual_positions[pos], this->actual_positions[pos + d]);
106 2486 sjacqu01
                        ++it_begin;
107 2486 sjacqu01
                        ++pos;
108 2486 sjacqu01
                    }
109 2486 sjacqu01
110 2486 sjacqu01
                    // calculate correct initial actual positions
111 2486 sjacqu01
                    for (std::size_t i = 1; i < 5; ++i)
112 2486 sjacqu01
                    {
113 2486 sjacqu01
                        this->actual_positions[i] += this->actual_positions[i - 1];
114 2486 sjacqu01
                    }
115 2486 sjacqu01
                }
116 2486 sjacqu01
            }
117 2486 sjacqu01
            else
118 2486 sjacqu01
            {
119 2486 sjacqu01
                std::size_t sample_cell = 1; // k
120 2486 sjacqu01
121 2486 sjacqu01
                // find cell k such that heights[k-1] <= args[sample] < heights[k] and adjust extreme values
122 2486 sjacqu01
                if (args[sample] < this->heights[0])
123 2486 sjacqu01
                {
124 2486 sjacqu01
                    this->heights[0] = args[sample];
125 2486 sjacqu01
                    this->actual_positions[0] = args[weight];
126 2486 sjacqu01
                    sample_cell = 1;
127 2486 sjacqu01
                }
128 2486 sjacqu01
                else if (this->heights[4] <= args[sample])
129 2486 sjacqu01
                {
130 2486 sjacqu01
                    this->heights[4] = args[sample];
131 2486 sjacqu01
                    sample_cell = 4;
132 2486 sjacqu01
                }
133 2486 sjacqu01
                else
134 2486 sjacqu01
                {
135 2486 sjacqu01
                    typedef typename array_type::iterator iterator;
136 2486 sjacqu01
                    iterator it = std::upper_bound(
137 2486 sjacqu01
                        this->heights.begin()
138 2486 sjacqu01
                      , this->heights.end()
139 2486 sjacqu01
                      , args[sample]
140 2486 sjacqu01
                    );
141 2486 sjacqu01
142 2486 sjacqu01
                    sample_cell = std::distance(this->heights.begin(), it);
143 2486 sjacqu01
                }
144 2486 sjacqu01
145 2486 sjacqu01
                // increment positions of markers above sample_cell
146 2486 sjacqu01
                for (std::size_t i = sample_cell; i < 5; ++i)
147 2486 sjacqu01
                {
148 2486 sjacqu01
                    this->actual_positions[i] += args[weight];
149 2486 sjacqu01
                }
150 2486 sjacqu01
151 2486 sjacqu01
                // update desired positions for all markers
152 2486 sjacqu01
                this->desired_positions[0] = this->actual_positions[0];
153 2486 sjacqu01
                this->desired_positions[1] = (sum_of_weights(args) - this->actual_positions[0])
154 2486 sjacqu01
                                           * this->p/2. + this->actual_positions[0];
155 2486 sjacqu01
                this->desired_positions[2] = (sum_of_weights(args) - this->actual_positions[0])
156 2486 sjacqu01
                                           * this->p + this->actual_positions[0];
157 2486 sjacqu01
                this->desired_positions[3] = (sum_of_weights(args) - this->actual_positions[0])
158 2486 sjacqu01
                                           * (1. + this->p)/2. + this->actual_positions[0];
159 2486 sjacqu01
                this->desired_positions[4] = sum_of_weights(args);
160 2486 sjacqu01
161 2486 sjacqu01
                // adjust height and actual positions of markers 1 to 3 if necessary
162 2486 sjacqu01
                for (std::size_t i = 1; i <= 3; ++i)
163 2486 sjacqu01
                {
164 2486 sjacqu01
                    // offset to desired positions
165 2486 sjacqu01
                    float_type d = this->desired_positions[i] - this->actual_positions[i];
166 2486 sjacqu01
167 2486 sjacqu01
                    // offset to next position
168 2486 sjacqu01
                    float_type dp = this->actual_positions[i + 1] - this->actual_positions[i];
169 2486 sjacqu01
170 2486 sjacqu01
                    // offset to previous position
171 2486 sjacqu01
                    float_type dm = this->actual_positions[i - 1] - this->actual_positions[i];
172 2486 sjacqu01
173 2486 sjacqu01
                    // height ds
174 2486 sjacqu01
                    float_type hp = (this->heights[i + 1] - this->heights[i]) / dp;
175 2486 sjacqu01
                    float_type hm = (this->heights[i - 1] - this->heights[i]) / dm;
176 2486 sjacqu01
177 2486 sjacqu01
                    if ( ( d >= 1. && dp > 1. ) || ( d <= -1. && dm < -1. ) )
178 2486 sjacqu01
                    {
179 2486 sjacqu01
                        short sign_d = static_cast<short>(d / std::abs(d));
180 2486 sjacqu01
181 2486 sjacqu01
                        // try adjusting heights[i] using p-squared formula
182 2486 sjacqu01
                        float_type h = this->heights[i] + sign_d / (dp - dm) * ( (sign_d - dm) * hp + (dp - sign_d) * hm );
183 2486 sjacqu01
184 2486 sjacqu01
                        if ( this->heights[i - 1] < h && h < this->heights[i + 1] )
185 2486 sjacqu01
                        {
186 2486 sjacqu01
                            this->heights[i] = h;
187 2486 sjacqu01
                        }
188 2486 sjacqu01
                        else
189 2486 sjacqu01
                        {
190 2486 sjacqu01
                            // use linear formula
191 2486 sjacqu01
                            if (d>0)
192 2486 sjacqu01
                            {
193 2486 sjacqu01
                                this->heights[i] += hp;
194 2486 sjacqu01
                            }
195 2486 sjacqu01
                            if (d<0)
196 2486 sjacqu01
                            {
197 2486 sjacqu01
                                this->heights[i] -= hm;
198 2486 sjacqu01
                            }
199 2486 sjacqu01
                        }
200 2486 sjacqu01
                        this->actual_positions[i] += sign_d;
201 2486 sjacqu01
                    }
202 2486 sjacqu01
                }
203 2486 sjacqu01
            }
204 2486 sjacqu01
        }
205 2486 sjacqu01
206 2486 sjacqu01
        result_type result(dont_care) const
207 2486 sjacqu01
        {
208 2486 sjacqu01
            return this->heights[2];
209 2486 sjacqu01
        }
210 2486 sjacqu01
211 2486 sjacqu01
    private:
212 2486 sjacqu01
        float_type p;                    // the quantile probability p
213 2486 sjacqu01
        array_type heights;              // q_i
214 2486 sjacqu01
        array_type actual_positions;     // n_i
215 2486 sjacqu01
        array_type desired_positions;    // n'_i
216 2486 sjacqu01
    };
217 2486 sjacqu01
218 2486 sjacqu01
} // namespace impl
219 2486 sjacqu01
220 2486 sjacqu01
///////////////////////////////////////////////////////////////////////////////
221 2486 sjacqu01
// tag::weighted_p_square_quantile
222 2486 sjacqu01
//
223 2486 sjacqu01
namespace tag
224 2486 sjacqu01
{
225 2486 sjacqu01
    struct weighted_p_square_quantile
226 2486 sjacqu01
      : depends_on<count, sum_of_weights>
227 2486 sjacqu01
    {
228 2486 sjacqu01
        typedef accumulators::impl::weighted_p_square_quantile_impl<mpl::_1, mpl::_2, regular> impl;
229 2486 sjacqu01
    };
230 2486 sjacqu01
    struct weighted_p_square_quantile_for_median
231 2486 sjacqu01
      : depends_on<count, sum_of_weights>
232 2486 sjacqu01
    {
233 2486 sjacqu01
        typedef accumulators::impl::weighted_p_square_quantile_impl<mpl::_1, mpl::_2, for_median> impl;
234 2486 sjacqu01
    };
235 2486 sjacqu01
}
236 2486 sjacqu01
237 2486 sjacqu01
///////////////////////////////////////////////////////////////////////////////
238 2486 sjacqu01
// extract::weighted_p_square_quantile
239 2486 sjacqu01
// extract::weighted_p_square_quantile_for_median
240 2486 sjacqu01
//
241 2486 sjacqu01
namespace extract
242 2486 sjacqu01
{
243 2486 sjacqu01
    extractor<tag::weighted_p_square_quantile> const weighted_p_square_quantile = {};
244 2486 sjacqu01
    extractor<tag::weighted_p_square_quantile_for_median> const weighted_p_square_quantile_for_median = {};
245 2486 sjacqu01
246 2486 sjacqu01
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_p_square_quantile)
247 2486 sjacqu01
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_p_square_quantile_for_median)
248 2486 sjacqu01
}
249 2486 sjacqu01
250 2486 sjacqu01
using extract::weighted_p_square_quantile;
251 2486 sjacqu01
using extract::weighted_p_square_quantile_for_median;
252 2486 sjacqu01
253 2486 sjacqu01
}} // namespace boost::accumulators
254 2486 sjacqu01
255 2486 sjacqu01
#endif