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

source code of boost/libs/accumulators/include/boost/accumulators/statistics/weighted_peaks_over_threshold.hpp