Statistics
| Revision:

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

History | View | Annotate | Download (11.9 kB)

1
///////////////////////////////////////////////////////////////////////////////
2
// extended_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_EXTENDED_SINGLE_QUANTILE_HPP_DE_01_01_2006
9
#define BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_QUANTILE_HPP_DE_01_01_2006
10

    
11
#include <vector>
12
#include <functional>
13
#include <boost/throw_exception.hpp>
14
#include <boost/range/begin.hpp>
15
#include <boost/range/end.hpp>
16
#include <boost/range/iterator_range.hpp>
17
#include <boost/iterator/transform_iterator.hpp>
18
#include <boost/iterator/counting_iterator.hpp>
19
#include <boost/iterator/permutation_iterator.hpp>
20
#include <boost/parameter/keyword.hpp>
21
#include <boost/mpl/placeholders.hpp>
22
#include <boost/type_traits/is_same.hpp>
23
#include <boost/accumulators/framework/accumulator_base.hpp>
24
#include <boost/accumulators/framework/extractor.hpp>
25
#include <boost/accumulators/numeric/functional.hpp>
26
#include <boost/accumulators/framework/parameters/sample.hpp>
27
#include <boost/accumulators/framework/depends_on.hpp>
28
#include <boost/accumulators/statistics_fwd.hpp>
29
#include <boost/accumulators/statistics/count.hpp>
30
#include <boost/accumulators/statistics/parameters/quantile_probability.hpp>
31
#include <boost/accumulators/statistics/extended_p_square.hpp>
32
#include <boost/accumulators/statistics/weighted_extended_p_square.hpp>
33
#include <boost/accumulators/statistics/times2_iterator.hpp>
34

    
35
#ifdef _MSC_VER
36
# pragma warning(push)
37
# pragma warning(disable: 4127) // conditional expression is constant
38
#endif
39

    
40
namespace boost { namespace accumulators
41
{
42

    
43
namespace impl
44
{
45
    ///////////////////////////////////////////////////////////////////////////////
46
    // extended_p_square_quantile_impl
47
    //  single quantile estimation
48
    /**
49
        @brief Quantile estimation using the extended \f$P^2\f$ algorithm for weighted and unweighted samples
50

51
        Uses the quantile estimates calculated by the extended \f$P^2\f$ algorithm to compute
52
        intermediate quantile estimates by means of quadratic interpolation.
53

54
        @param quantile_probability The probability of the quantile to be estimated.
55
    */
56
    template<typename Sample, typename Impl1, typename Impl2> // Impl1: weighted/unweighted // Impl2: linear/quadratic
57
    struct extended_p_square_quantile_impl
58
      : accumulator_base
59
    {
60
        typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
61
        typedef std::vector<float_type> array_type;
62
        typedef iterator_range<
63
            detail::lvalue_index_iterator<
64
                permutation_iterator<
65
                    typename array_type::const_iterator
66
                  , detail::times2_iterator
67
                >
68
            >
69
        > range_type;
70
        // for boost::result_of
71
        typedef float_type result_type;
72

    
73
        template<typename Args>
74
        extended_p_square_quantile_impl(Args const &args)
75
          : probabilities(
76
                boost::begin(args[extended_p_square_probabilities])
77
              , boost::end(args[extended_p_square_probabilities])
78
            )
79
        {
80
        }
81

    
82
        template<typename Args>
83
        result_type result(Args const &args) const
84
        {
85
            typedef
86
                typename mpl::if_<
87
                    is_same<Impl1, weighted>
88
                  , tag::weighted_extended_p_square
89
                  , tag::extended_p_square
90
                >::type
91
            extended_p_square_tag;
92

    
93
            extractor<extended_p_square_tag> const some_extended_p_square = {};
94

    
95
            array_type heights(some_extended_p_square(args).size());
96
            std::copy(some_extended_p_square(args).begin(), some_extended_p_square(args).end(), heights.begin());
97

    
98
            this->probability = args[quantile_probability];
99

    
100
            typename array_type::const_iterator iter_probs = std::lower_bound(this->probabilities.begin(), this->probabilities.end(), this->probability);
101
            std::size_t dist = std::distance(this->probabilities.begin(), iter_probs);
102
            typename array_type::const_iterator iter_heights = heights.begin() + dist;
103

    
104
            // If this->probability is not in a valid range return NaN or throw exception
105
            if (this->probability < *this->probabilities.begin() || this->probability > *(this->probabilities.end() - 1))
106
            {
107
                if (std::numeric_limits<result_type>::has_quiet_NaN)
108
                {
109
                    return std::numeric_limits<result_type>::quiet_NaN();
110
                }
111
                else
112
                {
113
                    std::ostringstream msg;
114
                    msg << "probability = " << this->probability << " is not in valid range (";
115
                    msg << *this->probabilities.begin() << ", " << *(this->probabilities.end() - 1) << ")";
116
                    boost::throw_exception(std::runtime_error(msg.str()));
117
                    return Sample(0);
118
                }
119

    
120
            }
121

    
122
            if (*iter_probs == this->probability)
123
            {
124
                return heights[dist];
125
            }
126
            else
127
            {
128
                result_type res;
129

    
130
                if (is_same<Impl2, linear>::value)
131
                {
132
                    /////////////////////////////////////////////////////////////////////////////////
133
                    // LINEAR INTERPOLATION
134
                    //
135
                    float_type p1 = *iter_probs;
136
                    float_type p0 = *(iter_probs - 1);
137
                    float_type h1 = *iter_heights;
138
                    float_type h0 = *(iter_heights - 1);
139

    
140
                    float_type a = numeric::fdiv(h1 - h0, p1 - p0);
141
                    float_type b = h1 - p1 * a;
142

    
143
                    res = a * this->probability + b;
144
                }
145
                else
146
                {
147
                    /////////////////////////////////////////////////////////////////////////////////
148
                    // QUADRATIC INTERPOLATION
149
                    //
150
                    float_type p0, p1, p2;
151
                    float_type h0, h1, h2;
152

    
153
                    if ( (dist == 1 || *iter_probs - this->probability <= this->probability - *(iter_probs - 1) ) && dist != this->probabilities.size() - 1 )
154
                    {
155
                        p0 = *(iter_probs - 1);
156
                        p1 = *iter_probs;
157
                        p2 = *(iter_probs + 1);
158
                        h0 = *(iter_heights - 1);
159
                        h1 = *iter_heights;
160
                        h2 = *(iter_heights + 1);
161
                    }
162
                    else
163
                    {
164
                        p0 = *(iter_probs - 2);
165
                        p1 = *(iter_probs - 1);
166
                        p2 = *iter_probs;
167
                        h0 = *(iter_heights - 2);
168
                        h1 = *(iter_heights - 1);
169
                        h2 = *iter_heights;
170
                    }
171

    
172
                    float_type hp21 = numeric::fdiv(h2 - h1, p2 - p1);
173
                    float_type hp10 = numeric::fdiv(h1 - h0, p1 - p0);
174
                    float_type p21  = numeric::fdiv(p2 * p2 - p1 * p1, p2 - p1);
175
                    float_type p10  = numeric::fdiv(p1 * p1 - p0 * p0, p1 - p0);
176

    
177
                    float_type a = numeric::fdiv(hp21 - hp10, p21 - p10);
178
                    float_type b = hp21 - a * p21;
179
                    float_type c = h2 - a * p2 * p2 - b * p2;
180

    
181
                    res = a * this->probability * this-> probability + b * this->probability + c;
182
                }
183

    
184
                return res;
185
            }
186

    
187
        }
188
    private:
189

    
190
        array_type probabilities;
191
        mutable float_type probability;
192

    
193
    };
194

    
195
} // namespace impl
196

    
197
///////////////////////////////////////////////////////////////////////////////
198
// tag::extended_p_square_quantile
199
//
200
namespace tag
201
{
202
    struct extended_p_square_quantile
203
      : depends_on<extended_p_square>
204
    {
205
        typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, unweighted, linear> impl;
206
    };
207
    struct extended_p_square_quantile_quadratic
208
      : depends_on<extended_p_square>
209
    {
210
        typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, unweighted, quadratic> impl;
211
    };
212
    struct weighted_extended_p_square_quantile
213
      : depends_on<weighted_extended_p_square>
214
    {
215
        typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, weighted, linear> impl;
216
    };
217
    struct weighted_extended_p_square_quantile_quadratic
218
      : depends_on<weighted_extended_p_square>
219
    {
220
        typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, weighted, quadratic> impl;
221
    };
222
}
223

    
224
///////////////////////////////////////////////////////////////////////////////
225
// extract::extended_p_square_quantile
226
// extract::weighted_extended_p_square_quantile
227
//
228
namespace extract
229
{
230
    extractor<tag::extended_p_square_quantile> const extended_p_square_quantile = {};
231
    extractor<tag::extended_p_square_quantile_quadratic> const extended_p_square_quantile_quadratic = {};
232
    extractor<tag::weighted_extended_p_square_quantile> const weighted_extended_p_square_quantile = {};
233
    extractor<tag::weighted_extended_p_square_quantile_quadratic> const weighted_extended_p_square_quantile_quadratic = {};
234

    
235
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_quantile)
236
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_quantile_quadratic)
237
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_extended_p_square_quantile)
238
    BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_extended_p_square_quantile_quadratic)
239
}
240

    
241
using extract::extended_p_square_quantile;
242
using extract::extended_p_square_quantile_quadratic;
243
using extract::weighted_extended_p_square_quantile;
244
using extract::weighted_extended_p_square_quantile_quadratic;
245

    
246
// extended_p_square_quantile(linear) -> extended_p_square_quantile
247
template<>
248
struct as_feature<tag::extended_p_square_quantile(linear)>
249
{
250
    typedef tag::extended_p_square_quantile type;
251
};
252

    
253
// extended_p_square_quantile(quadratic) -> extended_p_square_quantile_quadratic
254
template<>
255
struct as_feature<tag::extended_p_square_quantile(quadratic)>
256
{
257
    typedef tag::extended_p_square_quantile_quadratic type;
258
};
259

    
260
// weighted_extended_p_square_quantile(linear) -> weighted_extended_p_square_quantile
261
template<>
262
struct as_feature<tag::weighted_extended_p_square_quantile(linear)>
263
{
264
    typedef tag::weighted_extended_p_square_quantile type;
265
};
266

    
267
// weighted_extended_p_square_quantile(quadratic) -> weighted_extended_p_square_quantile_quadratic
268
template<>
269
struct as_feature<tag::weighted_extended_p_square_quantile(quadratic)>
270
{
271
    typedef tag::weighted_extended_p_square_quantile_quadratic type;
272
};
273

    
274
// for the purposes of feature-based dependency resolution,
275
// extended_p_square_quantile and weighted_extended_p_square_quantile
276
// provide the same feature as quantile
277
template<>
278
struct feature_of<tag::extended_p_square_quantile>
279
  : feature_of<tag::quantile>
280
{
281
};
282
template<>
283
struct feature_of<tag::extended_p_square_quantile_quadratic>
284
  : feature_of<tag::quantile>
285
{
286
};
287
// So that extended_p_square_quantile can be automatically substituted with
288
// weighted_extended_p_square_quantile when the weight parameter is non-void
289
template<>
290
struct as_weighted_feature<tag::extended_p_square_quantile>
291
{
292
    typedef tag::weighted_extended_p_square_quantile type;
293
};
294

    
295
template<>
296
struct feature_of<tag::weighted_extended_p_square_quantile>
297
  : feature_of<tag::extended_p_square_quantile>
298
{
299
};
300

    
301
// So that extended_p_square_quantile_quadratic can be automatically substituted with
302
// weighted_extended_p_square_quantile_quadratic when the weight parameter is non-void
303
template<>
304
struct as_weighted_feature<tag::extended_p_square_quantile_quadratic>
305
{
306
    typedef tag::weighted_extended_p_square_quantile_quadratic type;
307
};
308
template<>
309
struct feature_of<tag::weighted_extended_p_square_quantile_quadratic>
310
  : feature_of<tag::extended_p_square_quantile_quadratic>
311
{
312
};
313

    
314
}} // namespace boost::accumulators
315

    
316
#ifdef _MSC_VER
317
# pragma warning(pop)
318
#endif
319

    
320
#endif