accumulators/statistics/extended_p_square_quantile.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // extended_p_square_quantile.hpp
0003 //
0004 //  Copyright 2005 Daniel Egloff. Distributed under the Boost
0005 //  Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_QUANTILE_HPP_DE_01_01_2006
0009 #define BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_QUANTILE_HPP_DE_01_01_2006
0010
0011 #include <vector>
0012 #include <functional>
0013 #include <boost/throw_exception.hpp>
0014 #include <boost/range/begin.hpp>
0015 #include <boost/range/end.hpp>
0016 #include <boost/range/iterator_range.hpp>
0017 #include <boost/iterator/transform_iterator.hpp>
0018 #include <boost/iterator/counting_iterator.hpp>
0019 #include <boost/iterator/permutation_iterator.hpp>
0020 #include <boost/parameter/keyword.hpp>
0021 #include <boost/mpl/placeholders.hpp>
0022 #include <boost/type_traits/is_same.hpp>
0023 #include <boost/accumulators/framework/accumulator_base.hpp>
0024 #include <boost/accumulators/framework/extractor.hpp>
0025 #include <boost/accumulators/numeric/functional.hpp>
0026 #include <boost/accumulators/framework/parameters/sample.hpp>
0027 #include <boost/accumulators/framework/depends_on.hpp>
0028 #include <boost/accumulators/statistics_fwd.hpp>
0029 #include <boost/accumulators/statistics/count.hpp>
0030 #include <boost/accumulators/statistics/parameters/quantile_probability.hpp>
0031 #include <boost/accumulators/statistics/extended_p_square.hpp>
0032 #include <boost/accumulators/statistics/weighted_extended_p_square.hpp>
0033 #include <boost/accumulators/statistics/times2_iterator.hpp>
0034
0035 #ifdef _MSC_VER
0036 # pragma warning(push)
0037 # pragma warning(disable: 4127) // conditional expression is constant
0038 #endif
0039
0040 namespace boost { namespace accumulators
0041 {
0042
0043 namespace impl
0044 {
0045     ///////////////////////////////////////////////////////////////////////////////
0046     // extended_p_square_quantile_impl
0047     //  single quantile estimation
0048     /**
0049         @brief Quantile estimation using the extended \f$P^2\f$ algorithm for weighted and unweighted samples
0050
0051         Uses the quantile estimates calculated by the extended \f$P^2\f$ algorithm to compute
0052         intermediate quantile estimates by means of quadratic interpolation.
0053
0054         @param quantile_probability The probability of the quantile to be estimated.
0055     */
0056     template<typename Sample, typename Impl1, typename Impl2> // Impl1: weighted/unweighted // Impl2: linear/quadratic
0057     struct extended_p_square_quantile_impl
0058       : accumulator_base
0059     {
0060         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
0061         typedef std::vector<float_type> array_type;
0062         typedef iterator_range<
0063             detail::lvalue_index_iterator<
0064                 permutation_iterator<
0065                     typename array_type::const_iterator
0066                   , detail::times2_iterator
0067                 >
0068             >
0069         > range_type;
0070         // for boost::result_of
0071         typedef float_type result_type;
0072
0073         template<typename Args>
0074         extended_p_square_quantile_impl(Args const &args)
0075           : probabilities(
0076                 boost::begin(args[extended_p_square_probabilities])
0077               , boost::end(args[extended_p_square_probabilities])
0078             )
0079           , probability()
0080         {
0081         }
0082
0083         template<typename Args>
0084         result_type result(Args const &args) const
0085         {
0086             typedef
0087                 typename mpl::if_<
0088                     is_same<Impl1, weighted>
0089                   , tag::weighted_extended_p_square
0090                   , tag::extended_p_square
0091                 >::type
0092             extended_p_square_tag;
0093
0094             extractor<extended_p_square_tag> const some_extended_p_square = {};
0095
0096             array_type heights(some_extended_p_square(args).size());
0097             std::copy(some_extended_p_square(args).begin(), some_extended_p_square(args).end(), heights.begin());
0098
0099             this->probability = args[quantile_probability];
0100
0101             typename array_type::const_iterator iter_probs = std::lower_bound(this->probabilities.begin(), this->probabilities.end(), this->probability);
0102             std::size_t dist = std::distance(this->probabilities.begin(), iter_probs);
0103             typename array_type::const_iterator iter_heights = heights.begin() + dist;
0104
0105             // If this->probability is not in a valid range return NaN or throw exception
0106             if (this->probability < *this->probabilities.begin() || this->probability > *(this->probabilities.end() - 1))
0107             {
0108                 if (std::numeric_limits<result_type>::has_quiet_NaN)
0109                 {
0110                     return std::numeric_limits<result_type>::quiet_NaN();
0111                 }
0112                 else
0113                 {
0114                     std::ostringstream msg;
0115                     msg << "probability = " << this->probability << " is not in valid range (";
0116                     msg << *this->probabilities.begin() << ", " << *(this->probabilities.end() - 1) << ")";
0117                     boost::throw_exception(std::runtime_error(msg.str()));
0118                     return Sample(0);
0119                 }
0120
0121             }
0122
0123             if (*iter_probs == this->probability)
0124             {
0125                 return heights[dist];
0126             }
0127             else
0128             {
0129                 result_type res;
0130
0131                 if (is_same<Impl2, linear>::value)
0132                 {
0133                     /////////////////////////////////////////////////////////////////////////////////
0134                     // LINEAR INTERPOLATION
0135                     //
0136                     float_type p1 = *iter_probs;
0137                     float_type p0 = *(iter_probs - 1);
0138                     float_type h1 = *iter_heights;
0139                     float_type h0 = *(iter_heights - 1);
0140
0141                     float_type a = numeric::fdiv(h1 - h0, p1 - p0);
0142                     float_type b = h1 - p1 * a;
0143
0144                     res = a * this->probability + b;
0145                 }
0146                 else
0147                 {
0148                     /////////////////////////////////////////////////////////////////////////////////
0149                     // QUADRATIC INTERPOLATION
0150                     //
0151                     float_type p0, p1, p2;
0152                     float_type h0, h1, h2;
0153
0154                     if ( (dist == 1 || *iter_probs - this->probability <= this->probability - *(iter_probs - 1) ) && dist != this->probabilities.size() - 1 )
0155                     {
0156                         p0 = *(iter_probs - 1);
0157                         p1 = *iter_probs;
0158                         p2 = *(iter_probs + 1);
0159                         h0 = *(iter_heights - 1);
0160                         h1 = *iter_heights;
0161                         h2 = *(iter_heights + 1);
0162                     }
0163                     else
0164                     {
0165                         p0 = *(iter_probs - 2);
0166                         p1 = *(iter_probs - 1);
0167                         p2 = *iter_probs;
0168                         h0 = *(iter_heights - 2);
0169                         h1 = *(iter_heights - 1);
0170                         h2 = *iter_heights;
0171                     }
0172
0173                     float_type hp21 = numeric::fdiv(h2 - h1, p2 - p1);
0174                     float_type hp10 = numeric::fdiv(h1 - h0, p1 - p0);
0175                     float_type p21  = numeric::fdiv(p2 * p2 - p1 * p1, p2 - p1);
0176                     float_type p10  = numeric::fdiv(p1 * p1 - p0 * p0, p1 - p0);
0177
0178                     float_type a = numeric::fdiv(hp21 - hp10, p21 - p10);
0179                     float_type b = hp21 - a * p21;
0180                     float_type c = h2 - a * p2 * p2 - b * p2;
0181
0182                     res = a * this->probability * this-> probability + b * this->probability + c;
0183                 }
0184
0185                 return res;
0186             }
0187
0188         }
0189
0190     public:
0191         // make this accumulator serializeable
0192         // TODO: do we need to split to load/save and verify that the parameters did not change?
0193         template<class Archive>
0194         void serialize(Archive & ar, const unsigned int file_version)
0195         {
0196             ar & probabilities;
0197             ar & probability;
0198         }
0199
0200     private:
0201
0202         array_type probabilities;
0203         mutable float_type probability;
0204
0205     };
0206
0207 } // namespace impl
0208
0209 ///////////////////////////////////////////////////////////////////////////////
0210 // tag::extended_p_square_quantile
0211 //
0212 namespace tag
0213 {
0214     struct extended_p_square_quantile
0215       : depends_on<extended_p_square>
0216     {
0217         typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, unweighted, linear> impl;
0218     };
0219     struct extended_p_square_quantile_quadratic
0220       : depends_on<extended_p_square>
0221     {
0222         typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, unweighted, quadratic> impl;
0223     };
0224     struct weighted_extended_p_square_quantile
0225       : depends_on<weighted_extended_p_square>
0226     {
0227         typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, weighted, linear> impl;
0228     };
0229     struct weighted_extended_p_square_quantile_quadratic
0230       : depends_on<weighted_extended_p_square>
0231     {
0232         typedef accumulators::impl::extended_p_square_quantile_impl<mpl::_1, weighted, quadratic> impl;
0233     };
0234 }
0235
0236 ///////////////////////////////////////////////////////////////////////////////
0237 // extract::extended_p_square_quantile
0238 // extract::weighted_extended_p_square_quantile
0239 //
0240 namespace extract
0241 {
0242     extractor<tag::extended_p_square_quantile> const extended_p_square_quantile = {};
0243     extractor<tag::extended_p_square_quantile_quadratic> const extended_p_square_quantile_quadratic = {};
0244     extractor<tag::weighted_extended_p_square_quantile> const weighted_extended_p_square_quantile = {};
0245     extractor<tag::weighted_extended_p_square_quantile_quadratic> const weighted_extended_p_square_quantile_quadratic = {};
0246
0247     BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_quantile)
0248     BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square_quantile_quadratic)
0249     BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_extended_p_square_quantile)
0250     BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_extended_p_square_quantile_quadratic)
0251 }
0252
0253 using extract::extended_p_square_quantile;
0254 using extract::extended_p_square_quantile_quadratic;
0255 using extract::weighted_extended_p_square_quantile;
0256 using extract::weighted_extended_p_square_quantile_quadratic;
0257
0258 // extended_p_square_quantile(linear) -> extended_p_square_quantile
0259 template<>
0260 struct as_feature<tag::extended_p_square_quantile(linear)>
0261 {
0262     typedef tag::extended_p_square_quantile type;
0263 };
0264
0265 // extended_p_square_quantile(quadratic) -> extended_p_square_quantile_quadratic
0266 template<>
0267 struct as_feature<tag::extended_p_square_quantile(quadratic)>
0268 {
0269     typedef tag::extended_p_square_quantile_quadratic type;
0270 };
0271
0272 // weighted_extended_p_square_quantile(linear) -> weighted_extended_p_square_quantile
0273 template<>
0274 struct as_feature<tag::weighted_extended_p_square_quantile(linear)>
0275 {
0276     typedef tag::weighted_extended_p_square_quantile type;
0277 };
0278
0279 // weighted_extended_p_square_quantile(quadratic) -> weighted_extended_p_square_quantile_quadratic
0280 template<>
0281 struct as_feature<tag::weighted_extended_p_square_quantile(quadratic)>
0282 {
0283     typedef tag::weighted_extended_p_square_quantile_quadratic type;
0284 };
0285
0286 // for the purposes of feature-based dependency resolution,
0287 // extended_p_square_quantile and weighted_extended_p_square_quantile
0288 // provide the same feature as quantile
0289 template<>
0290 struct feature_of<tag::extended_p_square_quantile>
0291   : feature_of<tag::quantile>
0292 {
0293 };
0294 template<>
0295 struct feature_of<tag::extended_p_square_quantile_quadratic>
0296   : feature_of<tag::quantile>
0297 {
0298 };
0299 // So that extended_p_square_quantile can be automatically substituted with
0300 // weighted_extended_p_square_quantile when the weight parameter is non-void
0301 template<>
0302 struct as_weighted_feature<tag::extended_p_square_quantile>
0303 {
0304     typedef tag::weighted_extended_p_square_quantile type;
0305 };
0306
0307 template<>
0308 struct feature_of<tag::weighted_extended_p_square_quantile>
0309   : feature_of<tag::extended_p_square_quantile>
0310 {
0311 };
0312
0313 // So that extended_p_square_quantile_quadratic can be automatically substituted with
0314 // weighted_extended_p_square_quantile_quadratic when the weight parameter is non-void
0315 template<>
0316 struct as_weighted_feature<tag::extended_p_square_quantile_quadratic>
0317 {
0318     typedef tag::weighted_extended_p_square_quantile_quadratic type;
0319 };
0320 template<>
0321 struct feature_of<tag::weighted_extended_p_square_quantile_quadratic>
0322   : feature_of<tag::extended_p_square_quantile_quadratic>
0323 {
0324 };
0325
0326 }} // namespace boost::accumulators
0327
0328 #ifdef _MSC_VER
0329 # pragma warning(pop)
0330 #endif
0331
0332 #endif