/Src/Dependencies/Boost/boost/accumulators/statistics/extended_p_square.hpp
C++ Header | 293 lines | 195 code | 38 blank | 60 comment | 21 complexity | cd3bb37c768c007548bdfc8e8e30e6a6 MD5 | raw file
1/////////////////////////////////////////////////////////////////////////////// 2// extended_p_square.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_EXTENDED_SINGLE_HPP_DE_01_01_2006 9#define BOOST_ACCUMULATORS_STATISTICS_EXTENDED_SINGLE_HPP_DE_01_01_2006 10 11#include <vector> 12#include <functional> 13#include <boost/range/begin.hpp> 14#include <boost/range/end.hpp> 15#include <boost/range/iterator_range.hpp> 16#include <boost/iterator/transform_iterator.hpp> 17#include <boost/iterator/counting_iterator.hpp> 18#include <boost/iterator/permutation_iterator.hpp> 19#include <boost/parameter/keyword.hpp> 20#include <boost/mpl/placeholders.hpp> 21#include <boost/accumulators/framework/accumulator_base.hpp> 22#include <boost/accumulators/framework/extractor.hpp> 23#include <boost/accumulators/numeric/functional.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/count.hpp> 28#include <boost/accumulators/statistics/times2_iterator.hpp> 29 30namespace boost { namespace accumulators 31{ 32/////////////////////////////////////////////////////////////////////////////// 33// probabilities named parameter 34// 35BOOST_PARAMETER_NESTED_KEYWORD(tag, extended_p_square_probabilities, probabilities) 36 37namespace impl 38{ 39 /////////////////////////////////////////////////////////////////////////////// 40 // extended_p_square_impl 41 // multiple quantile estimation 42 /** 43 @brief Multiple quantile estimation with the extended \f$P^2\f$ algorithm 44 45 Extended \f$P^2\f$ algorithm for estimation of several quantiles without storing samples. 46 Assume that \f$m\f$ quantiles \f$\xi_{p_1}, \ldots, \xi_{p_m}\f$ are to be estimated. 47 Instead of storing the whole sample cumulative distribution, the algorithm maintains only 48 \f$m+2\f$ principal markers and \f$m+1\f$ middle markers, whose positions are updated 49 with each sample and whose heights are adjusted (if necessary) using a piecewise-parablic 50 formula. The heights of these central markers are the current estimates of the quantiles 51 and returned as an iterator range. 52 53 For further details, see 54 55 K. E. E. Raatikainen, Simultaneous estimation of several quantiles, Simulation, Volume 49, 56 Number 4 (October), 1986, p. 159-164. 57 58 The extended \f$ P^2 \f$ algorithm generalizess the \f$ P^2 \f$ algorithm of 59 60 R. Jain and I. Chlamtac, The P^2 algorithmus for dynamic calculation of quantiles and 61 histograms without storing observations, Communications of the ACM, 62 Volume 28 (October), Number 10, 1985, p. 1076-1085. 63 64 @param extended_p_square_probabilities A vector of quantile probabilities. 65 */ 66 template<typename Sample> 67 struct extended_p_square_impl 68 : accumulator_base 69 { 70 typedef typename numeric::functional::average<Sample, std::size_t>::result_type float_type; 71 typedef std::vector<float_type> array_type; 72 // for boost::result_of 73 typedef iterator_range< 74 detail::lvalue_index_iterator< 75 permutation_iterator< 76 typename array_type::const_iterator 77 , detail::times2_iterator 78 > 79 > 80 > result_type; 81 82 template<typename Args> 83 extended_p_square_impl(Args const &args) 84 : probabilities( 85 boost::begin(args[extended_p_square_probabilities]) 86 , boost::end(args[extended_p_square_probabilities]) 87 ) 88 , heights(2 * probabilities.size() + 3) 89 , actual_positions(heights.size()) 90 , desired_positions(heights.size()) 91 , positions_increments(heights.size()) 92 { 93 std::size_t num_quantiles = this->probabilities.size(); 94 std::size_t num_markers = this->heights.size(); 95 96 for(std::size_t i = 0; i < num_markers; ++i) 97 { 98 this->actual_positions[i] = i + 1; 99 } 100 101 this->positions_increments[0] = 0.; 102 this->positions_increments[num_markers - 1] = 1.; 103 104 for(std::size_t i = 0; i < num_quantiles; ++i) 105 { 106 this->positions_increments[2 * i + 2] = probabilities[i]; 107 } 108 109 for(std::size_t i = 0; i <= num_quantiles; ++i) 110 { 111 this->positions_increments[2 * i + 1] = 112 0.5 * (this->positions_increments[2 * i] + this->positions_increments[2 * i + 2]); 113 } 114 115 for(std::size_t i = 0; i < num_markers; ++i) 116 { 117 this->desired_positions[i] = 1. + 2. * (num_quantiles + 1.) * this->positions_increments[i]; 118 } 119 } 120 121 template<typename Args> 122 void operator ()(Args const &args) 123 { 124 std::size_t cnt = count(args); 125 126 // m+2 principal markers and m+1 middle markers 127 std::size_t num_markers = 2 * this->probabilities.size() + 3; 128 129 // first accumulate num_markers samples 130 if(cnt <= num_markers) 131 { 132 this->heights[cnt - 1] = args[sample]; 133 134 // complete the initialization of heights by sorting 135 if(cnt == num_markers) 136 { 137 std::sort(this->heights.begin(), this->heights.end()); 138 } 139 } 140 else 141 { 142 std::size_t sample_cell = 1; 143 144 // find cell k = sample_cell such that heights[k-1] <= sample < heights[k] 145 if(args[sample] < this->heights[0]) 146 { 147 this->heights[0] = args[sample]; 148 sample_cell = 1; 149 } 150 else if(args[sample] >= this->heights[num_markers - 1]) 151 { 152 this->heights[num_markers - 1] = args[sample]; 153 sample_cell = num_markers - 1; 154 } 155 else 156 { 157 typedef typename array_type::iterator iterator; 158 iterator it = std::upper_bound( 159 this->heights.begin() 160 , this->heights.end() 161 , args[sample] 162 ); 163 164 sample_cell = std::distance(this->heights.begin(), it); 165 } 166 167 // update actual positions of all markers above sample_cell index 168 for(std::size_t i = sample_cell; i < num_markers; ++i) 169 { 170 ++this->actual_positions[i]; 171 } 172 173 // update desired positions of all markers 174 for(std::size_t i = 0; i < num_markers; ++i) 175 { 176 this->desired_positions[i] += this->positions_increments[i]; 177 } 178 179 // adjust heights and actual positions of markers 1 to num_markers-2 if necessary 180 for(std::size_t i = 1; i <= num_markers - 2; ++i) 181 { 182 // offset to desired position 183 float_type d = this->desired_positions[i] - this->actual_positions[i]; 184 185 // offset to next position 186 float_type dp = this->actual_positions[i+1] - this->actual_positions[i]; 187 188 // offset to previous position 189 float_type dm = this->actual_positions[i-1] - this->actual_positions[i]; 190 191 // height ds 192 float_type hp = (this->heights[i+1] - this->heights[i]) / dp; 193 float_type hm = (this->heights[i-1] - this->heights[i]) / dm; 194 195 if((d >= 1 && dp > 1) || (d <= -1 && dm < -1)) 196 { 197 short sign_d = static_cast<short>(d / std::abs(d)); 198 199 float_type h = this->heights[i] + sign_d / (dp - dm) * ((sign_d - dm)*hp 200 + (dp - sign_d) * hm); 201 202 // try adjusting heights[i] using p-squared formula 203 if(this->heights[i - 1] < h && h < this->heights[i + 1]) 204 { 205 this->heights[i] = h; 206 } 207 else 208 { 209 // use linear formula 210 if(d > 0) 211 { 212 this->heights[i] += hp; 213 } 214 if(d < 0) 215 { 216 this->heights[i] -= hm; 217 } 218 } 219 this->actual_positions[i] += sign_d; 220 } 221 } 222 } 223 } 224 225 result_type result(dont_care) const 226 { 227 // for i in [1,probabilities.size()], return heights[i * 2] 228 detail::times2_iterator idx_begin = detail::make_times2_iterator(1); 229 detail::times2_iterator idx_end = detail::make_times2_iterator(this->probabilities.size() + 1); 230 231 return result_type( 232 make_permutation_iterator(this->heights.begin(), idx_begin) 233 , make_permutation_iterator(this->heights.begin(), idx_end) 234 ); 235 } 236 237 private: 238 array_type probabilities; // the quantile probabilities 239 array_type heights; // q_i 240 array_type actual_positions; // n_i 241 array_type desired_positions; // d_i 242 array_type positions_increments; // f_i 243 }; 244 245} // namespace impl 246 247/////////////////////////////////////////////////////////////////////////////// 248// tag::extended_p_square 249// 250namespace tag 251{ 252 struct extended_p_square 253 : depends_on<count> 254 , extended_p_square_probabilities 255 { 256 typedef accumulators::impl::extended_p_square_impl<mpl::_1> impl; 257 258 #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED 259 /// tag::extended_p_square::probabilities named paramter 260 static boost::parameter::keyword<tag::probabilities> const probabilities; 261 #endif 262 }; 263} 264 265/////////////////////////////////////////////////////////////////////////////// 266// extract::extended_p_square 267// 268namespace extract 269{ 270 extractor<tag::extended_p_square> const extended_p_square = {}; 271 272 BOOST_ACCUMULATORS_IGNORE_GLOBAL(extended_p_square) 273} 274 275using extract::extended_p_square; 276 277// So that extended_p_square can be automatically substituted with 278// weighted_extended_p_square when the weight parameter is non-void 279template<> 280struct as_weighted_feature<tag::extended_p_square> 281{ 282 typedef tag::weighted_extended_p_square type; 283}; 284 285template<> 286struct feature_of<tag::weighted_extended_p_square> 287 : feature_of<tag::extended_p_square> 288{ 289}; 290 291}} // namespace boost::accumulators 292 293#endif