boost/accumulators/statistics/weighted_peaks_over_threshold.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 // weighted_peaks_over_threshold.hpp
   3 //
   4 //  Copyright 2006 Daniel Egloff, Olivier Gygi. 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_WEIGHTED_PEAKS_OVER_THRESHOLD_HPP_DE_01_01_2006
   9 #define BOOST_ACCUMULATORS_STATISTICS_WEIGHTED_PEAKS_OVER_THRESHOLD_HPP_DE_01_01_2006
  10
  11 #include <vector>
  12 #include <limits>
  13 #include <numeric>
  14 #include <functional>
  15 #include <boost/range.hpp>
  16 #include <boost/mpl/if.hpp>
  17 #include <boost/mpl/placeholders.hpp>
  18 #include <boost/parameter/keyword.hpp>
  19 #include <boost/tuple/tuple.hpp>
  20 #include <boost/accumulators/numeric/functional.hpp>
  21 #include <boost/accumulators/framework/accumulator_base.hpp>
  22 #include <boost/accumulators/framework/extractor.hpp>
  23 #include <boost/accumulators/framework/parameters/sample.hpp>
  24 #include <boost/accumulators/framework/depends_on.hpp>
  25 #include <boost/accumulators/statistics_fwd.hpp>
  26 #include <boost/accumulators/statistics/parameters/quantile_probability.hpp>
  27 #include <boost/accumulators/statistics/peaks_over_threshold.hpp> // for named parameters pot_threshold_value and pot_threshold_probability
  28 #include <boost/accumulators/statistics/sum.hpp>
  29 #include <boost/accumulators/statistics/tail_variate.hpp>
  30
  31 #ifdef _MSC_VER
  32 # pragma warning(push)
  33 # pragma warning(disable: 4127) // conditional expression is constant
  34 #endif
  35
  36 namespace boost { namespace accumulators
  37 {
  38
  39 namespace impl
  40 {
  41
  42     ///////////////////////////////////////////////////////////////////////////////
  43     // weighted_peaks_over_threshold_impl
  44     //  works with an explicit threshold value and does not depend on order statistics of weighted samples
  45     /**
  46         @brief Weighted Peaks over Threshold Method for Weighted Quantile and Weighted Tail Mean Estimation
  47
  48         @sa peaks_over_threshold_impl
  49
  50         @param quantile_probability
  51         @param pot_threshold_value
  52     */
  53     template<typename Sample, typename Weight, typename LeftRight>
  54     struct weighted_peaks_over_threshold_impl
  55       : accumulator_base
  56     {
  57         typedef typename numeric::functional::multiplies<Weight, Sample>::result_type weighted_sample;
  58         typedef typename numeric::functional::average<weighted_sample, std::size_t>::result_type float_type;
  59         // for boost::result_of
  60         typedef boost::tuple<float_type, float_type, float_type> result_type;
  61
  62         template<typename Args>
  63         weighted_peaks_over_threshold_impl(Args const &args)
  64           : sign_((is_same<LeftRight, left>::value) ? -1 : 1)
  65           , mu_(sign_ * numeric::average(args[sample | Sample()], (std::size_t)1))
  66           , sigma2_(numeric::average(args[sample | Sample()], (std::size_t)1))
  67           , w_sum_(numeric::average(args[weight | Weight()], (std::size_t)1))
  68           , threshold_(sign_ * args[pot_threshold_value])
  69           , fit_parameters_(boost::make_tuple(0., 0., 0.))
  70           , is_dirty_(true)
  71         {
  72         }
  73
  74         template<typename Args>
  75         void operator ()(Args const &args)
  76         {
  77             this->is_dirty_ = true;
  78
  79             if (this->sign_ * args[sample] > this->threshold_)
  80             {
  81                 this->mu_ += args[weight] * args[sample];
  82                 this->sigma2_ += args[weight] * args[sample] * args[sample];
  83                 this->w_sum_ += args[weight];
  84             }
  85         }
  86
  87         template<typename Args>
  88         result_type result(Args const &args) const
  89         {
  90             if (this->is_dirty_)
  91             {
  92                 this->is_dirty_ = false;
  93
  94                 this->mu_ = this->sign_ * numeric::average(this->mu_, this->w_sum_);
  95                 this->sigma2_ = numeric::average(this->sigma2_, this->w_sum_);
  96                 this->sigma2_ -= this->mu_ * this->mu_;
  97
  98                 float_type threshold_probability = numeric::average(sum_of_weights(args) - this->w_sum_, sum_of_weights(args));
  99
 100                 float_type tmp = numeric::average(( this->mu_ - this->threshold_ )*( this->mu_ - this->threshold_ ), this->sigma2_);
 101                 float_type xi_hat = 0.5 * ( 1. - tmp );
 102                 float_type beta_hat = 0.5 * ( this->mu_ - this->threshold_ ) * ( 1. + tmp );
 103                 float_type beta_bar = beta_hat * std::pow(1. - threshold_probability, xi_hat);
 104                 float_type u_bar = this->threshold_ - beta_bar * ( std::pow(1. - threshold_probability, -xi_hat) - 1.)/xi_hat;
 105                 this->fit_parameters_ = boost::make_tuple(u_bar, beta_bar, xi_hat);
 106             }
 107
 108             return this->fit_parameters_;
 109         }
 110
 111     private:
 112         short sign_;                         // for left tail fitting, mirror the extreme values
 113         mutable float_type mu_;              // mean of samples above threshold
 114         mutable float_type sigma2_;          // variance of samples above threshold
 115         mutable float_type w_sum_;           // sum of weights of samples above threshold
 116         float_type threshold_;
 117         mutable result_type fit_parameters_; // boost::tuple that stores fit parameters
 118         mutable bool is_dirty_;
 119     };
 120
 121     ///////////////////////////////////////////////////////////////////////////////
 122     // weighted_peaks_over_threshold_prob_impl
 123     //  determines threshold from a given threshold probability using order statistics
 124     /**
 125         @brief Peaks over Threshold Method for Quantile and Tail Mean Estimation
 126
 127         @sa weighted_peaks_over_threshold_impl
 128
 129         @param quantile_probability
 130         @param pot_threshold_probability
 131     */
 132     template<typename Sample, typename Weight, typename LeftRight>
 133     struct weighted_peaks_over_threshold_prob_impl
 134       : accumulator_base
 135     {
 136         typedef typename numeric::functional::multiplies<Weight, Sample>::result_type weighted_sample;
 137         typedef typename numeric::functional::average<weighted_sample, std::size_t>::result_type float_type;
 138         // for boost::result_of
 139         typedef boost::tuple<float_type, float_type, float_type> result_type;
 140
 141         template<typename Args>
 142         weighted_peaks_over_threshold_prob_impl(Args const &args)
 143           : sign_((is_same<LeftRight, left>::value) ? -1 : 1)
 144           , mu_(sign_ * numeric::average(args[sample | Sample()], (std::size_t)1))
 145           , sigma2_(numeric::average(args[sample | Sample()], (std::size_t)1))
 146           , threshold_probability_(args[pot_threshold_probability])
 147           , fit_parameters_(boost::make_tuple(0., 0., 0.))
 148           , is_dirty_(true)
 149         {
 150         }
 151
 152         void operator ()(dont_care)
 153         {
 154             this->is_dirty_ = true;
 155         }
 156
 157         template<typename Args>
 158         result_type result(Args const &args) const
 159         {
 160             if (this->is_dirty_)
 161             {
 162                 this->is_dirty_ = false;
 163
 164                 float_type threshold = sum_of_weights(args)
 165                              * ( ( is_same<LeftRight, left>::value ) ? this->threshold_probability_ : 1. - this->threshold_probability_ );
 166
 167                 std::size_t n = 0;
 168                 Weight sum = Weight(0);
 169
 170                 while (sum < threshold)
 171                 {
 172                     if (n < static_cast<std::size_t>(tail_weights(args).size()))
 173                     {
 174                         mu_ += *(tail_weights(args).begin() + n) * *(tail(args).begin() + n);
 175                         sigma2_ += *(tail_weights(args).begin() + n) * *(tail(args).begin() + n) * (*(tail(args).begin() + n));
 176                         sum += *(tail_weights(args).begin() + n);
 177                         n++;
 178                     }
 179                     else
 180                     {
 181                         if (std::numeric_limits<float_type>::has_quiet_NaN)
 182                         {
 183                             return boost::make_tuple(
 184                                 std::numeric_limits<float_type>::quiet_NaN()
 185                               , std::numeric_limits<float_type>::quiet_NaN()
 186                               , std::numeric_limits<float_type>::quiet_NaN()
 187                             );
 188                         }
 189                         else
 190                         {
 191                             std::ostringstream msg;
 192                             msg << "index n = " << n << " is not in valid range [0, " << tail(args).size() << ")";
 193                             boost::throw_exception(std::runtime_error(msg.str()));
 194                             return boost::make_tuple(Sample(0), Sample(0), Sample(0));
 195                         }
 196                     }
 197                 }
 198
 199                 float_type u = *(tail(args).begin() + n - 1) * this->sign_;
 200
 201
 202                 this->mu_ = this->sign_ * numeric::average(this->mu_, sum);
 203                 this->sigma2_ = numeric::average(this->sigma2_, sum);
 204                 this->sigma2_ -= this->mu_ * this->mu_;
 205
 206                 if (is_same<LeftRight, left>::value)
 207                     this->threshold_probability_ = 1. - this->threshold_probability_;
 208
 209                 float_type tmp = numeric::average(( this->mu_ - u )*( this->mu_ - u ), this->sigma2_);
 210                 float_type xi_hat = 0.5 * ( 1. - tmp );
 211                 float_type beta_hat = 0.5 * ( this->mu_ - u ) * ( 1. + tmp );
 212                 float_type beta_bar = beta_hat * std::pow(1. - threshold_probability_, xi_hat);
 213                 float_type u_bar = u - beta_bar * ( std::pow(1. - threshold_probability_, -xi_hat) - 1.)/xi_hat;
 214                 this->fit_parameters_ = boost::make_tuple(u_bar, beta_bar, xi_hat);
 215
 216             }
 217
 218             return this->fit_parameters_;
 219         }
 220
 221     private:
 222         short sign_;                                // for left tail fitting, mirror the extreme values
 223         mutable float_type mu_;                     // mean of samples above threshold u
 224         mutable float_type sigma2_;                 // variance of samples above threshold u
 225         mutable float_type threshold_probability_;
 226         mutable result_type fit_parameters_;        // boost::tuple that stores fit parameters
 227         mutable bool is_dirty_;
 228     };
 229
 230 } // namespace impl
 231
 232 ///////////////////////////////////////////////////////////////////////////////
 233 // tag::weighted_peaks_over_threshold
 234 //
 235 namespace tag
 236 {
 237     template<typename LeftRight>
 238     struct weighted_peaks_over_threshold
 239       : depends_on<sum_of_weights>
 240       , pot_threshold_value
 241     {
 242         /// INTERNAL ONLY
 243         typedef accumulators::impl::weighted_peaks_over_threshold_impl<mpl::_1, mpl::_2, LeftRight> impl;
 244     };
 245
 246     template<typename LeftRight>
 247     struct weighted_peaks_over_threshold_prob
 248       : depends_on<sum_of_weights, tail_weights<LeftRight> >
 249       , pot_threshold_probability
 250     {
 251         /// INTERNAL ONLY
 252         typedef accumulators::impl::weighted_peaks_over_threshold_prob_impl<mpl::_1, mpl::_2, LeftRight> impl;
 253     };
 254 }
 255
 256 ///////////////////////////////////////////////////////////////////////////////
 257 // extract::weighted_peaks_over_threshold
 258 //
 259 namespace extract
 260 {
 261     extractor<tag::abstract_peaks_over_threshold> const weighted_peaks_over_threshold = {};
 262 }
 263
 264 using extract::weighted_peaks_over_threshold;
 265
 266 // weighted_peaks_over_threshold<LeftRight>(with_threshold_value) -> weighted_peaks_over_threshold<LeftRight>
 267 template<typename LeftRight>
 268 struct as_feature<tag::weighted_peaks_over_threshold<LeftRight>(with_threshold_value)>
 269 {
 270     typedef tag::weighted_peaks_over_threshold<LeftRight> type;
 271 };
 272
 273 // weighted_peaks_over_threshold<LeftRight>(with_threshold_probability) -> weighted_peaks_over_threshold_prob<LeftRight>
 274 template<typename LeftRight>
 275 struct as_feature<tag::weighted_peaks_over_threshold<LeftRight>(with_threshold_probability)>
 276 {
 277     typedef tag::weighted_peaks_over_threshold_prob<LeftRight> type;
 278 };
 279
 280 }} // namespace boost::accumulators
 281
 282 #ifdef _MSC_VER
 283 # pragma warning(pop)
 284 #endif
 285
 286 #endif