boost/accumulators/statistics/variance.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 // 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
   8 #ifndef BOOST_ACCUMULATORS_STATISTICS_VARIANCE_HPP_EAN_28_10_2005
   9 #define BOOST_ACCUMULATORS_STATISTICS_VARIANCE_HPP_EAN_28_10_2005
  10
  11 #include <boost/mpl/placeholders.hpp>
  12 #include <boost/accumulators/framework/accumulator_base.hpp>
  13 #include <boost/accumulators/framework/extractor.hpp>
  14 #include <boost/accumulators/numeric/functional.hpp>
  15 #include <boost/accumulators/framework/parameters/sample.hpp>
  16 #include <boost/accumulators/framework/depends_on.hpp>
  17 #include <boost/accumulators/statistics_fwd.hpp>
  18 #include <boost/accumulators/statistics/count.hpp>
  19 #include <boost/accumulators/statistics/sum.hpp>
  20 #include <boost/accumulators/statistics/mean.hpp>
  21 #include <boost/accumulators/statistics/moment.hpp>
  22
  23 namespace boost { namespace accumulators
  24 {
  25
  26 namespace impl
  27 {
  28     //! Lazy calculation of variance.
  29     /*!
  30         Default sample variance implementation based on the second moment \f$ M_n^{(2)} \f$ moment<2>, mean and count.
  31         \f[
  32             \sigma_n^2 = M_n^{(2)} - \mu_n^2.
  33         \f]
  34         where
  35         \f[
  36             \mu_n = \frac{1}{n} \sum_{i = 1}^n x_i.
  37         \f]
  38         is the estimate of the sample mean and \f$n\f$ is the number of samples.
  39     */
  40     template<typename Sample, typename MeanFeature>
  41     struct lazy_variance_impl
  42       : accumulator_base
  43     {
  44         // for boost::result_of
  45         typedef typename numeric::functional::average<Sample, std::size_t>::result_type result_type;
  46
  47         lazy_variance_impl(dont_care) {}
  48
  49         template<typename Args>
  50         result_type result(Args const &args) const
  51         {
  52             extractor<MeanFeature> mean;
  53             result_type tmp = mean(args);
  54             return moment<2>(args) - tmp * tmp;
  55         }
  56     };
  57
  58     //! Iterative calculation of variance.
  59     /*!
  60         Iterative calculation of sample variance \f$\sigma_n^2\f$ according to the formula
  61         \f[
  62             \sigma_n^2 = \frac{1}{n} \sum_{i = 1}^n (x_i - \mu_n)^2 = \frac{n-1}{n} \sigma_{n-1}^2 + \frac{1}{n-1}(x_n - \mu_n)^2.
  63         \f]
  64         where
  65         \f[
  66             \mu_n = \frac{1}{n} \sum_{i = 1}^n x_i.
  67         \f]
  68         is the estimate of the sample mean and \f$n\f$ is the number of samples.
  69
  70         Note that the sample variance is not defined for \f$n <= 1\f$.
  71
  72         A simplification can be obtained by the approximate recursion
  73         \f[
  74             \sigma_n^2 \approx \frac{n-1}{n} \sigma_{n-1}^2 + \frac{1}{n}(x_n - \mu_n)^2.
  75         \f]
  76         because the difference
  77         \f[
  78             \left(\frac{1}{n-1} - \frac{1}{n}\right)(x_n - \mu_n)^2 = \frac{1}{n(n-1)}(x_n - \mu_n)^2.
  79         \f]
  80         converges to zero as \f$n \rightarrow \infty\f$. However, for small \f$ n \f$ the difference
  81         can be non-negligible.
  82     */
  83     template<typename Sample, typename MeanFeature, typename Tag>
  84     struct variance_impl
  85       : accumulator_base
  86     {
  87         // for boost::result_of
  88         typedef typename numeric::functional::average<Sample, std::size_t>::result_type result_type;
  89
  90         template<typename Args>
  91         variance_impl(Args const &args)
  92           : variance(numeric::average(args[sample | Sample()], numeric::one<std::size_t>::value))
  93         {
  94         }
  95
  96         template<typename Args>
  97         void operator ()(Args const &args)
  98         {
  99             std::size_t cnt = count(args);
 100
 101             if(cnt > 1)
 102             {
 103                 extractor<MeanFeature> mean;
 104                 result_type tmp = args[parameter::keyword<Tag>::get()] - mean(args);
 105                 this->variance =
 106                     numeric::average(this->variance * (cnt - 1), cnt)
 107                   + numeric::average(tmp * tmp, cnt - 1);
 108             }
 109         }
 110
 111         result_type result(dont_care) const
 112         {
 113             return this->variance;
 114         }
 115
 116     private:
 117         result_type variance;
 118     };
 119
 120 } // namespace impl
 121
 122 ///////////////////////////////////////////////////////////////////////////////
 123 // tag::variance
 124 // tag::immediate_variance
 125 //
 126 namespace tag
 127 {
 128     struct lazy_variance
 129       : depends_on<moment<2>, mean>
 130     {
 131         /// INTERNAL ONLY
 132         ///
 133         typedef accumulators::impl::lazy_variance_impl<mpl::_1, mean> impl;
 134     };
 135
 136     struct variance
 137       : depends_on<count, immediate_mean>
 138     {
 139         /// INTERNAL ONLY
 140         ///
 141         typedef accumulators::impl::variance_impl<mpl::_1, mean, sample> impl;
 142     };
 143 }
 144
 145 ///////////////////////////////////////////////////////////////////////////////
 146 // extract::lazy_variance
 147 // extract::variance
 148 //
 149 namespace extract
 150 {
 151     extractor<tag::lazy_variance> const lazy_variance = {};
 152     extractor<tag::variance> const variance = {};
 153 }
 154
 155 using extract::lazy_variance;
 156 using extract::variance;
 157
 158 // variance(lazy) -> lazy_variance
 159 template<>
 160 struct as_feature<tag::variance(lazy)>
 161 {
 162     typedef tag::lazy_variance type;
 163 };
 164
 165 // variance(immediate) -> variance
 166 template<>
 167 struct as_feature<tag::variance(immediate)>
 168 {
 169     typedef tag::variance type;
 170 };
 171
 172 // for the purposes of feature-based dependency resolution,
 173 // immediate_variance provides the same feature as variance
 174 template<>
 175 struct feature_of<tag::lazy_variance>
 176   : feature_of<tag::variance>
 177 {
 178 };
 179
 180 // So that variance can be automatically substituted with
 181 // weighted_variance when the weight parameter is non-void.
 182 template<>
 183 struct as_weighted_feature<tag::variance>
 184 {
 185     typedef tag::weighted_variance type;
 186 };
 187
 188 // for the purposes of feature-based dependency resolution,
 189 // weighted_variance provides the same feature as variance
 190 template<>
 191 struct feature_of<tag::weighted_variance>
 192   : feature_of<tag::variance>
 193 {
 194 };
 195
 196 // So that immediate_variance can be automatically substituted with
 197 // immediate_weighted_variance when the weight parameter is non-void.
 198 template<>
 199 struct as_weighted_feature<tag::lazy_variance>
 200 {
 201     typedef tag::lazy_weighted_variance type;
 202 };
 203
 204 // for the purposes of feature-based dependency resolution,
 205 // immediate_weighted_variance provides the same feature as immediate_variance
 206 template<>
 207 struct feature_of<tag::lazy_weighted_variance>
 208   : feature_of<tag::lazy_variance>
 209 {
 210 };
 211
 212 ////////////////////////////////////////////////////////////////////////////
 213 //// droppable_accumulator<variance_impl>
 214 ////  need to specialize droppable lazy variance to cache the result at the
 215 ////  point the accumulator is dropped.
 216 ///// INTERNAL ONLY
 217 /////
 218 //template<typename Sample, typename MeanFeature>
 219 //struct droppable_accumulator<impl::variance_impl<Sample, MeanFeature> >
 220 //  : droppable_accumulator_base<
 221 //        with_cached_result<impl::variance_impl<Sample, MeanFeature> >
 222 //    >
 223 //{
 224 //    template<typename Args>
 225 //    droppable_accumulator(Args const &args)
 226 //      : droppable_accumulator::base(args)
 227 //    {
 228 //    }
 229 //};
 230
 231 }} // namespace boost::accumulators
 232
 233 #endif