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

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