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

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