/Src/Dependencies/Boost/boost/accumulators/statistics/weighted_peaks_over_threshold.hpp

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