accumulators/statistics/median.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // median.hpp
0003 //
0004 //  Copyright 2006 Eric Niebler, Olivier Gygi. 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_MEDIAN_HPP_EAN_28_10_2005
0009 #define BOOST_ACCUMULATORS_STATISTICS_MEDIAN_HPP_EAN_28_10_2005
0010
0011 #include <boost/mpl/placeholders.hpp>
0012 #include <boost/range/iterator_range.hpp>
0013 #include <boost/accumulators/framework/accumulator_base.hpp>
0014 #include <boost/accumulators/framework/extractor.hpp>
0015 #include <boost/accumulators/numeric/functional.hpp>
0016 #include <boost/accumulators/framework/parameters/sample.hpp>
0017 #include <boost/accumulators/framework/depends_on.hpp>
0018 #include <boost/accumulators/statistics_fwd.hpp>
0019 #include <boost/accumulators/statistics/count.hpp>
0020 #include <boost/accumulators/statistics/p_square_quantile.hpp>
0021 #include <boost/accumulators/statistics/density.hpp>
0022 #include <boost/accumulators/statistics/p_square_cumul_dist.hpp>
0023
0024 namespace boost { namespace accumulators
0025 {
0026
0027 namespace impl
0028 {
0029     ///////////////////////////////////////////////////////////////////////////////
0030     // median_impl
0031     //
0032     /**
0033         @brief Median estimation based on the \f$P^2\f$ quantile estimator
0034
0035         The \f$P^2\f$ algorithm is invoked with a quantile probability of 0.5.
0036     */
0037     template<typename Sample>
0038     struct median_impl
0039       : accumulator_base
0040     {
0041         // for boost::result_of
0042         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type result_type;
0043
0044         median_impl(dont_care) {}
0045
0046         template<typename Args>
0047         result_type result(Args const &args) const
0048         {
0049             return p_square_quantile_for_median(args);
0050         }
0051
0052         // serialization is done by accumulators it depends on
0053         template<class Archive>
0054         void serialize(Archive & ar, const unsigned int file_version) {}
0055     };
0056     ///////////////////////////////////////////////////////////////////////////////
0057     // with_density_median_impl
0058     //
0059     /**
0060         @brief Median estimation based on the density estimator
0061
0062         The algorithm determines the bin in which the \f$0.5*cnt\f$-th sample lies, \f$cnt\f$ being
0063         the total number of samples. It returns the approximate horizontal position of this sample,
0064         based on a linear interpolation inside the bin.
0065     */
0066     template<typename Sample>
0067     struct with_density_median_impl
0068       : accumulator_base
0069     {
0070         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
0071         typedef std::vector<std::pair<float_type, float_type> > histogram_type;
0072         typedef iterator_range<typename histogram_type::iterator> range_type;
0073         // for boost::result_of
0074         typedef float_type result_type;
0075
0076         template<typename Args>
0077         with_density_median_impl(Args const &args)
0078           : sum(numeric::fdiv(args[sample | Sample()], (std::size_t)1))
0079           , is_dirty(true)
0080         {
0081         }
0082
0083         void operator ()(dont_care)
0084         {
0085             this->is_dirty = true;
0086         }
0087
0088
0089         template<typename Args>
0090         result_type result(Args const &args) const
0091         {
0092             if (this->is_dirty)
0093             {
0094                 this->is_dirty = false;
0095
0096                 std::size_t cnt = count(args);
0097                 range_type histogram = density(args);
0098                 typename range_type::iterator it = histogram.begin();
0099                 while (this->sum < 0.5 * cnt)
0100                 {
0101                     this->sum += it->second * cnt;
0102                     ++it;
0103                 }
0104                 --it;
0105                 float_type over = numeric::fdiv(this->sum - 0.5 * cnt, it->second * cnt);
0106                 this->median = it->first * over + (it + 1)->first * (1. - over);
0107             }
0108
0109             return this->median;
0110         }
0111
0112         // make this accumulator serializeable
0113         template<class Archive>
0114         void serialize(Archive & ar, const unsigned int file_version)
0115         {
0116             ar & sum;
0117             ar & is_dirty;
0118             ar & median;
0119         }
0120
0121
0122     private:
0123         mutable float_type sum;
0124         mutable bool is_dirty;
0125         mutable float_type median;
0126     };
0127
0128     ///////////////////////////////////////////////////////////////////////////////
0129     // with_p_square_cumulative_distribution_median_impl
0130     //
0131     /**
0132         @brief Median estimation based on the \f$P^2\f$ cumulative distribution estimator
0133
0134         The algorithm determines the first (leftmost) bin with a height exceeding 0.5. It
0135         returns the approximate horizontal position of where the cumulative distribution
0136         equals 0.5, based on a linear interpolation inside the bin.
0137     */
0138     template<typename Sample>
0139     struct with_p_square_cumulative_distribution_median_impl
0140       : accumulator_base
0141     {
0142         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type float_type;
0143         typedef std::vector<std::pair<float_type, float_type> > histogram_type;
0144         typedef iterator_range<typename histogram_type::iterator> range_type;
0145         // for boost::result_of
0146         typedef float_type result_type;
0147
0148         with_p_square_cumulative_distribution_median_impl(dont_care)
0149           : is_dirty(true)
0150         {
0151         }
0152
0153         void operator ()(dont_care)
0154         {
0155             this->is_dirty = true;
0156         }
0157
0158         template<typename Args>
0159         result_type result(Args const &args) const
0160         {
0161             if (this->is_dirty)
0162             {
0163                 this->is_dirty = false;
0164
0165                 range_type histogram = p_square_cumulative_distribution(args);
0166                 typename range_type::iterator it = histogram.begin();
0167                 while (it->second < 0.5)
0168                 {
0169                     ++it;
0170                 }
0171                 float_type over = numeric::fdiv(it->second - 0.5, it->second - (it - 1)->second);
0172                 this->median = it->first * over + (it + 1)->first * ( 1. - over );
0173             }
0174
0175             return this->median;
0176         }
0177
0178         // make this accumulator serializeable
0179         template<class Archive>
0180         void serialize(Archive & ar, const unsigned int file_version)
0181         {
0182             ar & is_dirty;
0183             ar & median;
0184         }
0185
0186     private:
0187
0188         mutable bool is_dirty;
0189         mutable float_type median;
0190     };
0191
0192 } // namespace impl
0193
0194 ///////////////////////////////////////////////////////////////////////////////
0195 // tag::median
0196 // tag::with_densisty_median
0197 // tag::with_p_square_cumulative_distribution_median
0198 //
0199 namespace tag
0200 {
0201     struct median
0202       : depends_on<p_square_quantile_for_median>
0203     {
0204         /// INTERNAL ONLY
0205         ///
0206         typedef accumulators::impl::median_impl<mpl::_1> impl;
0207     };
0208     struct with_density_median
0209       : depends_on<count, density>
0210     {
0211         /// INTERNAL ONLY
0212         ///
0213         typedef accumulators::impl::with_density_median_impl<mpl::_1> impl;
0214     };
0215     struct with_p_square_cumulative_distribution_median
0216       : depends_on<p_square_cumulative_distribution>
0217     {
0218         /// INTERNAL ONLY
0219         ///
0220         typedef accumulators::impl::with_p_square_cumulative_distribution_median_impl<mpl::_1> impl;
0221     };
0222 }
0223
0224 ///////////////////////////////////////////////////////////////////////////////
0225 // extract::median
0226 // extract::with_density_median
0227 // extract::with_p_square_cumulative_distribution_median
0228 //
0229 namespace extract
0230 {
0231     extractor<tag::median> const median = {};
0232     extractor<tag::with_density_median> const with_density_median = {};
0233     extractor<tag::with_p_square_cumulative_distribution_median> const with_p_square_cumulative_distribution_median = {};
0234
0235     BOOST_ACCUMULATORS_IGNORE_GLOBAL(median)
0236     BOOST_ACCUMULATORS_IGNORE_GLOBAL(with_density_median)
0237     BOOST_ACCUMULATORS_IGNORE_GLOBAL(with_p_square_cumulative_distribution_median)
0238 }
0239
0240 using extract::median;
0241 using extract::with_density_median;
0242 using extract::with_p_square_cumulative_distribution_median;
0243
0244 // median(with_p_square_quantile) -> median
0245 template<>
0246 struct as_feature<tag::median(with_p_square_quantile)>
0247 {
0248     typedef tag::median type;
0249 };
0250
0251 // median(with_density) -> with_density_median
0252 template<>
0253 struct as_feature<tag::median(with_density)>
0254 {
0255     typedef tag::with_density_median type;
0256 };
0257
0258 // median(with_p_square_cumulative_distribution) -> with_p_square_cumulative_distribution_median
0259 template<>
0260 struct as_feature<tag::median(with_p_square_cumulative_distribution)>
0261 {
0262     typedef tag::with_p_square_cumulative_distribution_median type;
0263 };
0264
0265 // for the purposes of feature-based dependency resolution,
0266 // with_density_median and with_p_square_cumulative_distribution_median
0267 // provide the same feature as median
0268 template<>
0269 struct feature_of<tag::with_density_median>
0270   : feature_of<tag::median>
0271 {
0272 };
0273
0274 template<>
0275 struct feature_of<tag::with_p_square_cumulative_distribution_median>
0276   : feature_of<tag::median>
0277 {
0278 };
0279
0280 // So that median can be automatically substituted with
0281 // weighted_median when the weight parameter is non-void.
0282 template<>
0283 struct as_weighted_feature<tag::median>
0284 {
0285     typedef tag::weighted_median type;
0286 };
0287
0288 template<>
0289 struct feature_of<tag::weighted_median>
0290   : feature_of<tag::median>
0291 {
0292 };
0293
0294 // So that with_density_median can be automatically substituted with
0295 // with_density_weighted_median when the weight parameter is non-void.
0296 template<>
0297 struct as_weighted_feature<tag::with_density_median>
0298 {
0299     typedef tag::with_density_weighted_median type;
0300 };
0301
0302 template<>
0303 struct feature_of<tag::with_density_weighted_median>
0304   : feature_of<tag::with_density_median>
0305 {
0306 };
0307
0308 // So that with_p_square_cumulative_distribution_median can be automatically substituted with
0309 // with_p_square_cumulative_distribution_weighted_median when the weight parameter is non-void.
0310 template<>
0311 struct as_weighted_feature<tag::with_p_square_cumulative_distribution_median>
0312 {
0313     typedef tag::with_p_square_cumulative_distribution_weighted_median type;
0314 };
0315
0316 template<>
0317 struct feature_of<tag::with_p_square_cumulative_distribution_weighted_median>
0318   : feature_of<tag::with_p_square_cumulative_distribution_median>
0319 {
0320 };
0321
0322 }} // namespace boost::accumulators
0323
0324 #endif