Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // variance.hpp
Chris@16: //
Chris@16: //  Copyright 2005 Daniel Egloff, Eric Niebler. Distributed under the Boost
Chris@16: //  Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_ACCUMULATORS_STATISTICS_VARIANCE_HPP_EAN_28_10_2005
Chris@16: #define BOOST_ACCUMULATORS_STATISTICS_VARIANCE_HPP_EAN_28_10_2005
Chris@16: 
Chris@16: #include <boost/mpl/placeholders.hpp>
Chris@16: #include <boost/accumulators/framework/accumulator_base.hpp>
Chris@16: #include <boost/accumulators/framework/extractor.hpp>
Chris@16: #include <boost/accumulators/numeric/functional.hpp>
Chris@16: #include <boost/accumulators/framework/parameters/sample.hpp>
Chris@16: #include <boost/accumulators/framework/depends_on.hpp>
Chris@16: #include <boost/accumulators/statistics_fwd.hpp>
Chris@16: #include <boost/accumulators/statistics/count.hpp>
Chris@16: #include <boost/accumulators/statistics/sum.hpp>
Chris@16: #include <boost/accumulators/statistics/mean.hpp>
Chris@16: #include <boost/accumulators/statistics/moment.hpp>
Chris@16: 
Chris@16: namespace boost { namespace accumulators
Chris@16: {
Chris@16: 
Chris@16: namespace impl
Chris@16: {
Chris@16:     //! Lazy calculation of variance.
Chris@16:     /*!
Chris@16:         Default sample variance implementation based on the second moment \f$ M_n^{(2)} \f$ moment<2>, mean and count.
Chris@16:         \f[
Chris@16:             \sigma_n^2 = M_n^{(2)} - \mu_n^2.
Chris@16:         \f]
Chris@16:         where
Chris@16:         \f[
Chris@16:             \mu_n = \frac{1}{n} \sum_{i = 1}^n x_i.
Chris@16:         \f]
Chris@16:         is the estimate of the sample mean and \f$n\f$ is the number of samples.
Chris@16:     */
Chris@16:     template<typename Sample, typename MeanFeature>
Chris@16:     struct lazy_variance_impl
Chris@16:       : accumulator_base
Chris@16:     {
Chris@16:         // for boost::result_of
Chris@16:         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type result_type;
Chris@16: 
Chris@16:         lazy_variance_impl(dont_care) {}
Chris@16: 
Chris@16:         template<typename Args>
Chris@16:         result_type result(Args const &args) const
Chris@16:         {
Chris@16:             extractor<MeanFeature> mean;
Chris@16:             result_type tmp = mean(args);
Chris@16:             return accumulators::moment<2>(args) - tmp * tmp;
Chris@16:         }
Chris@16:     };
Chris@16: 
Chris@16:     //! Iterative calculation of variance.
Chris@16:     /*!
Chris@16:         Iterative calculation of sample variance \f$\sigma_n^2\f$ according to the formula
Chris@16:         \f[
Chris@16:             \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.
Chris@16:         \f]
Chris@16:         where
Chris@16:         \f[
Chris@16:             \mu_n = \frac{1}{n} \sum_{i = 1}^n x_i.
Chris@16:         \f]
Chris@16:         is the estimate of the sample mean and \f$n\f$ is the number of samples.
Chris@16: 
Chris@16:         Note that the sample variance is not defined for \f$n <= 1\f$.
Chris@16: 
Chris@16:         A simplification can be obtained by the approximate recursion
Chris@16:         \f[
Chris@16:             \sigma_n^2 \approx \frac{n-1}{n} \sigma_{n-1}^2 + \frac{1}{n}(x_n - \mu_n)^2.
Chris@16:         \f]
Chris@16:         because the difference
Chris@16:         \f[
Chris@16:             \left(\frac{1}{n-1} - \frac{1}{n}\right)(x_n - \mu_n)^2 = \frac{1}{n(n-1)}(x_n - \mu_n)^2.
Chris@16:         \f]
Chris@16:         converges to zero as \f$n \rightarrow \infty\f$. However, for small \f$ n \f$ the difference
Chris@16:         can be non-negligible.
Chris@16:     */
Chris@16:     template<typename Sample, typename MeanFeature, typename Tag>
Chris@16:     struct variance_impl
Chris@16:       : accumulator_base
Chris@16:     {
Chris@16:         // for boost::result_of
Chris@16:         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type result_type;
Chris@16: 
Chris@16:         template<typename Args>
Chris@16:         variance_impl(Args const &args)
Chris@16:           : variance(numeric::fdiv(args[sample | Sample()], numeric::one<std::size_t>::value))
Chris@16:         {
Chris@16:         }
Chris@16: 
Chris@16:         template<typename Args>
Chris@16:         void operator ()(Args const &args)
Chris@16:         {
Chris@16:             std::size_t cnt = count(args);
Chris@16: 
Chris@16:             if(cnt > 1)
Chris@16:             {
Chris@16:                 extractor<MeanFeature> mean;
Chris@16:                 result_type tmp = args[parameter::keyword<Tag>::get()] - mean(args);
Chris@16:                 this->variance =
Chris@16:                     numeric::fdiv(this->variance * (cnt - 1), cnt)
Chris@16:                   + numeric::fdiv(tmp * tmp, cnt - 1);
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         result_type result(dont_care) const
Chris@16:         {
Chris@16:             return this->variance;
Chris@16:         }
Chris@16: 
Chris@16:     private:
Chris@16:         result_type variance;
Chris@16:     };
Chris@16: 
Chris@16: } // namespace impl
Chris@16: 
Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // tag::variance
Chris@16: // tag::immediate_variance
Chris@16: //
Chris@16: namespace tag
Chris@16: {
Chris@16:     struct lazy_variance
Chris@16:       : depends_on<moment<2>, mean>
Chris@16:     {
Chris@16:         /// INTERNAL ONLY
Chris@16:         ///
Chris@16:         typedef accumulators::impl::lazy_variance_impl<mpl::_1, mean> impl;
Chris@16:     };
Chris@16: 
Chris@16:     struct variance
Chris@16:       : depends_on<count, immediate_mean>
Chris@16:     {
Chris@16:         /// INTERNAL ONLY
Chris@16:         ///
Chris@16:         typedef accumulators::impl::variance_impl<mpl::_1, mean, sample> impl;
Chris@16:     };
Chris@16: }
Chris@16: 
Chris@16: ///////////////////////////////////////////////////////////////////////////////
Chris@16: // extract::lazy_variance
Chris@16: // extract::variance
Chris@16: //
Chris@16: namespace extract
Chris@16: {
Chris@16:     extractor<tag::lazy_variance> const lazy_variance = {};
Chris@16:     extractor<tag::variance> const variance = {};
Chris@16: 
Chris@16:     BOOST_ACCUMULATORS_IGNORE_GLOBAL(lazy_variance)
Chris@16:     BOOST_ACCUMULATORS_IGNORE_GLOBAL(variance)
Chris@16: }
Chris@16: 
Chris@16: using extract::lazy_variance;
Chris@16: using extract::variance;
Chris@16: 
Chris@16: // variance(lazy) -> lazy_variance
Chris@16: template<>
Chris@16: struct as_feature<tag::variance(lazy)>
Chris@16: {
Chris@16:     typedef tag::lazy_variance type;
Chris@16: };
Chris@16: 
Chris@16: // variance(immediate) -> variance
Chris@16: template<>
Chris@16: struct as_feature<tag::variance(immediate)>
Chris@16: {
Chris@16:     typedef tag::variance type;
Chris@16: };
Chris@16: 
Chris@16: // for the purposes of feature-based dependency resolution,
Chris@16: // immediate_variance provides the same feature as variance
Chris@16: template<>
Chris@16: struct feature_of<tag::lazy_variance>
Chris@16:   : feature_of<tag::variance>
Chris@16: {
Chris@16: };
Chris@16: 
Chris@16: // So that variance can be automatically substituted with
Chris@16: // weighted_variance when the weight parameter is non-void.
Chris@16: template<>
Chris@16: struct as_weighted_feature<tag::variance>
Chris@16: {
Chris@16:     typedef tag::weighted_variance type;
Chris@16: };
Chris@16: 
Chris@16: // for the purposes of feature-based dependency resolution,
Chris@16: // weighted_variance provides the same feature as variance
Chris@16: template<>
Chris@16: struct feature_of<tag::weighted_variance>
Chris@16:   : feature_of<tag::variance>
Chris@16: {
Chris@16: };
Chris@16: 
Chris@16: // So that immediate_variance can be automatically substituted with
Chris@16: // immediate_weighted_variance when the weight parameter is non-void.
Chris@16: template<>
Chris@16: struct as_weighted_feature<tag::lazy_variance>
Chris@16: {
Chris@16:     typedef tag::lazy_weighted_variance type;
Chris@16: };
Chris@16: 
Chris@16: // for the purposes of feature-based dependency resolution,
Chris@16: // immediate_weighted_variance provides the same feature as immediate_variance
Chris@16: template<>
Chris@16: struct feature_of<tag::lazy_weighted_variance>
Chris@16:   : feature_of<tag::lazy_variance>
Chris@16: {
Chris@16: };
Chris@16: 
Chris@16: ////////////////////////////////////////////////////////////////////////////
Chris@16: //// droppable_accumulator<variance_impl>
Chris@16: ////  need to specialize droppable lazy variance to cache the result at the
Chris@16: ////  point the accumulator is dropped.
Chris@16: ///// INTERNAL ONLY
Chris@16: /////
Chris@16: //template<typename Sample, typename MeanFeature>
Chris@16: //struct droppable_accumulator<impl::variance_impl<Sample, MeanFeature> >
Chris@16: //  : droppable_accumulator_base<
Chris@16: //        with_cached_result<impl::variance_impl<Sample, MeanFeature> >
Chris@16: //    >
Chris@16: //{
Chris@16: //    template<typename Args>
Chris@16: //    droppable_accumulator(Args const &args)
Chris@16: //      : droppable_accumulator::base(args)
Chris@16: //    {
Chris@16: //    }
Chris@16: //};
Chris@16: 
Chris@16: }} // namespace boost::accumulators
Chris@16: 
Chris@16: #endif