weighted_variance.hpp 6.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // weighted_variance.hpp
  3. //
  4. // Copyright 2005 Daniel Egloff, Eric Niebler. 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_VARIANCE_HPP_EAN_28_10_2005
  8. #define BOOST_ACCUMULATORS_STATISTICS_WEIGHTED_VARIANCE_HPP_EAN_28_10_2005
  9. #include <boost/mpl/placeholders.hpp>
  10. #include <boost/accumulators/framework/accumulator_base.hpp>
  11. #include <boost/accumulators/framework/extractor.hpp>
  12. #include <boost/accumulators/numeric/functional.hpp>
  13. #include <boost/accumulators/framework/parameters/sample.hpp>
  14. #include <boost/accumulators/framework/depends_on.hpp>
  15. #include <boost/accumulators/statistics_fwd.hpp>
  16. #include <boost/accumulators/statistics/count.hpp>
  17. #include <boost/accumulators/statistics/variance.hpp>
  18. #include <boost/accumulators/statistics/weighted_sum.hpp>
  19. #include <boost/accumulators/statistics/weighted_mean.hpp>
  20. #include <boost/accumulators/statistics/weighted_moment.hpp>
  21. namespace boost { namespace accumulators
  22. {
  23. namespace impl
  24. {
  25. //! Lazy calculation of variance of weighted samples.
  26. /*!
  27. The default implementation of the variance of weighted samples is based on the second moment
  28. \f$\widehat{m}_n^{(2)}\f$ (weighted_moment<2>) and the mean\f$ \hat{\mu}_n\f$ (weighted_mean):
  29. \f[
  30. \hat{\sigma}_n^2 = \widehat{m}_n^{(2)}-\hat{\mu}_n^2,
  31. \f]
  32. where \f$n\f$ is the number of samples.
  33. */
  34. template<typename Sample, typename Weight, typename MeanFeature>
  35. struct lazy_weighted_variance_impl
  36. : accumulator_base
  37. {
  38. typedef typename numeric::functional::multiplies<Sample, Weight>::result_type weighted_sample;
  39. // for boost::result_of
  40. typedef typename numeric::functional::fdiv<weighted_sample, Weight>::result_type result_type;
  41. lazy_weighted_variance_impl(dont_care) {}
  42. template<typename Args>
  43. result_type result(Args const &args) const
  44. {
  45. extractor<MeanFeature> const some_mean = {};
  46. result_type tmp = some_mean(args);
  47. return accumulators::weighted_moment<2>(args) - tmp * tmp;
  48. }
  49. };
  50. //! Iterative calculation of variance of weighted samples.
  51. /*!
  52. Iterative calculation of variance of weighted samples:
  53. \f[
  54. \hat{\sigma}_n^2 =
  55. \frac{\bar{w}_n - w_n}{\bar{w}_n}\hat{\sigma}_{n - 1}^2
  56. + \frac{w_n}{\bar{w}_n - w_n}\left(X_n - \hat{\mu}_n\right)^2
  57. ,\quad n\ge2,\quad\hat{\sigma}_0^2 = 0.
  58. \f]
  59. where \f$\bar{w}_n\f$ is the sum of the \f$n\f$ weights \f$w_i\f$ and \f$\hat{\mu}_n\f$
  60. the estimate of the mean of the weighted samples. Note that the sample variance is not defined for
  61. \f$n <= 1\f$.
  62. */
  63. template<typename Sample, typename Weight, typename MeanFeature, typename Tag>
  64. struct weighted_variance_impl
  65. : accumulator_base
  66. {
  67. typedef typename numeric::functional::multiplies<Sample, Weight>::result_type weighted_sample;
  68. // for boost::result_of
  69. typedef typename numeric::functional::fdiv<weighted_sample, Weight>::result_type result_type;
  70. template<typename Args>
  71. weighted_variance_impl(Args const &args)
  72. : weighted_variance(numeric::fdiv(args[sample | Sample()], numeric::one<Weight>::value))
  73. {
  74. }
  75. template<typename Args>
  76. void operator ()(Args const &args)
  77. {
  78. std::size_t cnt = count(args);
  79. if(cnt > 1)
  80. {
  81. extractor<MeanFeature> const some_mean = {};
  82. result_type tmp = args[parameter::keyword<Tag>::get()] - some_mean(args);
  83. this->weighted_variance =
  84. numeric::fdiv(this->weighted_variance * (sum_of_weights(args) - args[weight]), sum_of_weights(args))
  85. + numeric::fdiv(tmp * tmp * args[weight], sum_of_weights(args) - args[weight] );
  86. }
  87. }
  88. result_type result(dont_care) const
  89. {
  90. return this->weighted_variance;
  91. }
  92. // make this accumulator serializeable
  93. template<class Archive>
  94. void serialize(Archive & ar, const unsigned int file_version)
  95. {
  96. ar & weighted_variance;
  97. }
  98. private:
  99. result_type weighted_variance;
  100. };
  101. } // namespace impl
  102. ///////////////////////////////////////////////////////////////////////////////
  103. // tag::weighted_variance
  104. // tag::immediate_weighted_variance
  105. //
  106. namespace tag
  107. {
  108. struct lazy_weighted_variance
  109. : depends_on<weighted_moment<2>, weighted_mean>
  110. {
  111. /// INTERNAL ONLY
  112. ///
  113. typedef accumulators::impl::lazy_weighted_variance_impl<mpl::_1, mpl::_2, weighted_mean> impl;
  114. };
  115. struct weighted_variance
  116. : depends_on<count, immediate_weighted_mean>
  117. {
  118. /// INTERNAL ONLY
  119. ///
  120. typedef accumulators::impl::weighted_variance_impl<mpl::_1, mpl::_2, immediate_weighted_mean, sample> impl;
  121. };
  122. }
  123. ///////////////////////////////////////////////////////////////////////////////
  124. // extract::weighted_variance
  125. // extract::immediate_weighted_variance
  126. //
  127. namespace extract
  128. {
  129. extractor<tag::lazy_weighted_variance> const lazy_weighted_variance = {};
  130. extractor<tag::weighted_variance> const weighted_variance = {};
  131. BOOST_ACCUMULATORS_IGNORE_GLOBAL(lazy_weighted_variance)
  132. BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_variance)
  133. }
  134. using extract::lazy_weighted_variance;
  135. using extract::weighted_variance;
  136. // weighted_variance(lazy) -> lazy_weighted_variance
  137. template<>
  138. struct as_feature<tag::weighted_variance(lazy)>
  139. {
  140. typedef tag::lazy_weighted_variance type;
  141. };
  142. // weighted_variance(immediate) -> weighted_variance
  143. template<>
  144. struct as_feature<tag::weighted_variance(immediate)>
  145. {
  146. typedef tag::weighted_variance type;
  147. };
  148. ////////////////////////////////////////////////////////////////////////////
  149. //// droppable_accumulator<weighted_variance_impl>
  150. //// need to specialize droppable lazy weighted_variance to cache the result at the
  151. //// point the accumulator is dropped.
  152. ///// INTERNAL ONLY
  153. /////
  154. //template<typename Sample, typename Weight, typename MeanFeature>
  155. //struct droppable_accumulator<impl::weighted_variance_impl<Sample, Weight, MeanFeature> >
  156. // : droppable_accumulator_base<
  157. // with_cached_result<impl::weighted_variance_impl<Sample, Weight, MeanFeature> >
  158. // >
  159. //{
  160. // template<typename Args>
  161. // droppable_accumulator(Args const &args)
  162. // : droppable_accumulator::base(args)
  163. // {
  164. // }
  165. //};
  166. }} // namespace boost::accumulators
  167. #endif