/Src/Dependencies/Boost/boost/accumulators/statistics/weighted_peaks_over_threshold.hpp
C++ Header | 288 lines | 207 code | 41 blank | 40 comment | 7 complexity | 6e798679fe737043cf39db0849cc80c4 MD5 | raw file
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 36namespace boost { namespace accumulators 37{ 38 39namespace 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// 235namespace 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// 259namespace extract 260{ 261 extractor<tag::abstract_peaks_over_threshold> const weighted_peaks_over_threshold = {}; 262 263 BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_peaks_over_threshold) 264} 265 266using extract::weighted_peaks_over_threshold; 267 268// weighted_peaks_over_threshold<LeftRight>(with_threshold_value) -> weighted_peaks_over_threshold<LeftRight> 269template<typename LeftRight> 270struct as_feature<tag::weighted_peaks_over_threshold<LeftRight>(with_threshold_value)> 271{ 272 typedef tag::weighted_peaks_over_threshold<LeftRight> type; 273}; 274 275// weighted_peaks_over_threshold<LeftRight>(with_threshold_probability) -> weighted_peaks_over_threshold_prob<LeftRight> 276template<typename LeftRight> 277struct as_feature<tag::weighted_peaks_over_threshold<LeftRight>(with_threshold_probability)> 278{ 279 typedef tag::weighted_peaks_over_threshold_prob<LeftRight> type; 280}; 281 282}} // namespace boost::accumulators 283 284#ifdef _MSC_VER 285# pragma warning(pop) 286#endif 287 288#endif