/Src/Dependencies/Boost/boost/accumulators/statistics/weighted_p_square_quantile.hpp
C++ Header | 254 lines | 165 code | 33 blank | 56 comment | 19 complexity | ed80122eb899d127440ca636519a5315 MD5 | raw file
1/////////////////////////////////////////////////////////////////////////////// 2// weighted_p_square_quantile.hpp 3// 4// Copyright 2005 Daniel Egloff. 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_P_SQUARE_QUANTILE_HPP_DE_01_01_2006 9#define BOOST_ACCUMULATORS_STATISTICS_WEIGHTED_P_SQUARE_QUANTILE_HPP_DE_01_01_2006 10 11#include <functional> 12#include <boost/array.hpp> 13#include <boost/parameter/keyword.hpp> 14#include <boost/mpl/placeholders.hpp> 15#include <boost/type_traits/is_same.hpp> 16#include <boost/accumulators/framework/accumulator_base.hpp> 17#include <boost/accumulators/framework/extractor.hpp> 18#include <boost/accumulators/numeric/functional.hpp> 19#include <boost/accumulators/framework/parameters/sample.hpp> 20#include <boost/accumulators/statistics_fwd.hpp> 21#include <boost/accumulators/statistics/count.hpp> 22#include <boost/accumulators/statistics/sum.hpp> 23#include <boost/accumulators/statistics/parameters/quantile_probability.hpp> 24 25namespace boost { namespace accumulators 26{ 27 28namespace impl { 29 /////////////////////////////////////////////////////////////////////////////// 30 // weighted_p_square_quantile_impl 31 // single quantile estimation with weighted samples 32 /** 33 @brief Single quantile estimation with the \f$P^2\f$ algorithm for weighted samples 34 35 This version of the \f$P^2\f$ algorithm extends the \f$P^2\f$ algorithm to support weighted samples. 36 The \f$P^2\f$ algorithm estimates a quantile dynamically without storing samples. Instead of 37 storing the whole sample cumulative distribution, only five points (markers) are stored. The heights 38 of these markers are the minimum and the maximum of the samples and the current estimates of the 39 \f$(p/2)\f$-, \f$p\f$ - and \f$(1+p)/2\f$ -quantiles. Their positions are equal to the number 40 of samples that are smaller or equal to the markers. Each time a new sample is added, the 41 positions of the markers are updated and if necessary their heights are adjusted using a piecewise- 42 parabolic formula. 43 44 For further details, see 45 46 R. Jain and I. Chlamtac, The P^2 algorithmus for dynamic calculation of quantiles and 47 histograms without storing observations, Communications of the ACM, 48 Volume 28 (October), Number 10, 1985, p. 1076-1085. 49 50 @param quantile_probability 51 */ 52 template<typename Sample, typename Weight, typename Impl> 53 struct weighted_p_square_quantile_impl 54 : accumulator_base 55 { 56 typedef typename numeric::functional::multiplies<Sample, Weight>::result_type weighted_sample; 57 typedef typename numeric::functional::average<weighted_sample, std::size_t>::result_type float_type; 58 typedef array<float_type, 5> array_type; 59 // for boost::result_of 60 typedef float_type result_type; 61 62 template<typename Args> 63 weighted_p_square_quantile_impl(Args const &args) 64 : p(is_same<Impl, for_median>::value ? 0.5 : args[quantile_probability | 0.5]) 65 , heights() 66 , actual_positions() 67 , desired_positions() 68 { 69 } 70 71 template<typename Args> 72 void operator ()(Args const &args) 73 { 74 std::size_t cnt = count(args); 75 76 // accumulate 5 first samples 77 if (cnt <= 5) 78 { 79 this->heights[cnt - 1] = args[sample]; 80 81 // In this initialization phase, actual_positions stores the weights of the 82 // inital samples that are needed at the end of the initialization phase to 83 // compute the correct initial positions of the markers. 84 this->actual_positions[cnt - 1] = args[weight]; 85 86 // complete the initialization of heights and actual_positions by sorting 87 if (cnt == 5) 88 { 89 // TODO: we need to sort the initial samples (in heights) in ascending order and 90 // sort their weights (in actual_positions) the same way. The following lines do 91 // it, but there must be a better and more efficient way of doing this. 92 typename array_type::iterator it_begin, it_end, it_min; 93 94 it_begin = this->heights.begin(); 95 it_end = this->heights.end(); 96 97 std::size_t pos = 0; 98 99 while (it_begin != it_end) 100 { 101 it_min = std::min_element(it_begin, it_end); 102 std::size_t d = std::distance(it_begin, it_min); 103 std::swap(*it_begin, *it_min); 104 std::swap(this->actual_positions[pos], this->actual_positions[pos + d]); 105 ++it_begin; 106 ++pos; 107 } 108 109 // calculate correct initial actual positions 110 for (std::size_t i = 1; i < 5; ++i) 111 { 112 this->actual_positions[i] += this->actual_positions[i - 1]; 113 } 114 } 115 } 116 else 117 { 118 std::size_t sample_cell = 1; // k 119 120 // find cell k such that heights[k-1] <= args[sample] < heights[k] and adjust extreme values 121 if (args[sample] < this->heights[0]) 122 { 123 this->heights[0] = args[sample]; 124 this->actual_positions[0] = args[weight]; 125 sample_cell = 1; 126 } 127 else if (this->heights[4] <= args[sample]) 128 { 129 this->heights[4] = args[sample]; 130 sample_cell = 4; 131 } 132 else 133 { 134 typedef typename array_type::iterator iterator; 135 iterator it = std::upper_bound( 136 this->heights.begin() 137 , this->heights.end() 138 , args[sample] 139 ); 140 141 sample_cell = std::distance(this->heights.begin(), it); 142 } 143 144 // increment positions of markers above sample_cell 145 for (std::size_t i = sample_cell; i < 5; ++i) 146 { 147 this->actual_positions[i] += args[weight]; 148 } 149 150 // update desired positions for all markers 151 this->desired_positions[0] = this->actual_positions[0]; 152 this->desired_positions[1] = (sum_of_weights(args) - this->actual_positions[0]) 153 * this->p/2. + this->actual_positions[0]; 154 this->desired_positions[2] = (sum_of_weights(args) - this->actual_positions[0]) 155 * this->p + this->actual_positions[0]; 156 this->desired_positions[3] = (sum_of_weights(args) - this->actual_positions[0]) 157 * (1. + this->p)/2. + this->actual_positions[0]; 158 this->desired_positions[4] = sum_of_weights(args); 159 160 // adjust height and actual positions of markers 1 to 3 if necessary 161 for (std::size_t i = 1; i <= 3; ++i) 162 { 163 // offset to desired positions 164 float_type d = this->desired_positions[i] - this->actual_positions[i]; 165 166 // offset to next position 167 float_type dp = this->actual_positions[i + 1] - this->actual_positions[i]; 168 169 // offset to previous position 170 float_type dm = this->actual_positions[i - 1] - this->actual_positions[i]; 171 172 // height ds 173 float_type hp = (this->heights[i + 1] - this->heights[i]) / dp; 174 float_type hm = (this->heights[i - 1] - this->heights[i]) / dm; 175 176 if ( ( d >= 1. && dp > 1. ) || ( d <= -1. && dm < -1. ) ) 177 { 178 short sign_d = static_cast<short>(d / std::abs(d)); 179 180 // try adjusting heights[i] using p-squared formula 181 float_type h = this->heights[i] + sign_d / (dp - dm) * ( (sign_d - dm) * hp + (dp - sign_d) * hm ); 182 183 if ( this->heights[i - 1] < h && h < this->heights[i + 1] ) 184 { 185 this->heights[i] = h; 186 } 187 else 188 { 189 // use linear formula 190 if (d>0) 191 { 192 this->heights[i] += hp; 193 } 194 if (d<0) 195 { 196 this->heights[i] -= hm; 197 } 198 } 199 this->actual_positions[i] += sign_d; 200 } 201 } 202 } 203 } 204 205 result_type result(dont_care) const 206 { 207 return this->heights[2]; 208 } 209 210 private: 211 float_type p; // the quantile probability p 212 array_type heights; // q_i 213 array_type actual_positions; // n_i 214 array_type desired_positions; // n'_i 215 }; 216 217} // namespace impl 218 219/////////////////////////////////////////////////////////////////////////////// 220// tag::weighted_p_square_quantile 221// 222namespace tag 223{ 224 struct weighted_p_square_quantile 225 : depends_on<count, sum_of_weights> 226 { 227 typedef accumulators::impl::weighted_p_square_quantile_impl<mpl::_1, mpl::_2, regular> impl; 228 }; 229 struct weighted_p_square_quantile_for_median 230 : depends_on<count, sum_of_weights> 231 { 232 typedef accumulators::impl::weighted_p_square_quantile_impl<mpl::_1, mpl::_2, for_median> impl; 233 }; 234} 235 236/////////////////////////////////////////////////////////////////////////////// 237// extract::weighted_p_square_quantile 238// extract::weighted_p_square_quantile_for_median 239// 240namespace extract 241{ 242 extractor<tag::weighted_p_square_quantile> const weighted_p_square_quantile = {}; 243 extractor<tag::weighted_p_square_quantile_for_median> const weighted_p_square_quantile_for_median = {}; 244 245 BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_p_square_quantile) 246 BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_p_square_quantile_for_median) 247} 248 249using extract::weighted_p_square_quantile; 250using extract::weighted_p_square_quantile_for_median; 251 252}} // namespace boost::accumulators 253 254#endif