accumulators/statistics/variance.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // variance.hpp
0003 //
0004 //  Copyright 2005 Daniel Egloff, Eric Niebler. Distributed under the Boost
0005 //  Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_ACCUMULATORS_STATISTICS_VARIANCE_HPP_EAN_28_10_2005
0009 #define BOOST_ACCUMULATORS_STATISTICS_VARIANCE_HPP_EAN_28_10_2005
0010
0011 #include <boost/mpl/placeholders.hpp>
0012 #include <boost/accumulators/framework/accumulator_base.hpp>
0013 #include <boost/accumulators/framework/extractor.hpp>
0014 #include <boost/accumulators/numeric/functional.hpp>
0015 #include <boost/accumulators/framework/parameters/sample.hpp>
0016 #include <boost/accumulators/framework/depends_on.hpp>
0017 #include <boost/accumulators/statistics_fwd.hpp>
0018 #include <boost/accumulators/statistics/count.hpp>
0019 #include <boost/accumulators/statistics/sum.hpp>
0020 #include <boost/accumulators/statistics/mean.hpp>
0021 #include <boost/accumulators/statistics/moment.hpp>
0022
0023 namespace boost { namespace accumulators
0024 {
0025
0026 namespace impl
0027 {
0028     //! Lazy calculation of variance.
0029     /*!
0030         Default sample variance implementation based on the second moment \f$ M_n^{(2)} \f$ moment<2>, mean and count.
0031         \f[
0032             \sigma_n^2 = M_n^{(2)} - \mu_n^2.
0033         \f]
0034         where
0035         \f[
0036             \mu_n = \frac{1}{n} \sum_{i = 1}^n x_i.
0037         \f]
0038         is the estimate of the sample mean and \f$n\f$ is the number of samples.
0039     */
0040     template<typename Sample, typename MeanFeature>
0041     struct lazy_variance_impl
0042       : accumulator_base
0043     {
0044         // for boost::result_of
0045         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type result_type;
0046
0047         lazy_variance_impl(dont_care) {}
0048
0049         template<typename Args>
0050         result_type result(Args const &args) const
0051         {
0052             extractor<MeanFeature> mean;
0053             result_type tmp = mean(args);
0054             return accumulators::moment<2>(args) - tmp * tmp;
0055         }
0056
0057         // serialization is done by accumulators it depends on
0058         template<class Archive>
0059         void serialize(Archive & /* ar */, const unsigned int /* file_version */) {}
0060     };
0061
0062     //! Iterative calculation of variance.
0063     /*!
0064         Iterative calculation of sample variance \f$\sigma_n^2\f$ according to the formula
0065         \f[
0066             \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.
0067         \f]
0068         where
0069         \f[
0070             \mu_n = \frac{1}{n} \sum_{i = 1}^n x_i.
0071         \f]
0072         is the estimate of the sample mean and \f$n\f$ is the number of samples.
0073
0074         Note that the sample variance is not defined for \f$n <= 1\f$.
0075
0076         A simplification can be obtained by the approximate recursion
0077         \f[
0078             \sigma_n^2 \approx \frac{n-1}{n} \sigma_{n-1}^2 + \frac{1}{n}(x_n - \mu_n)^2.
0079         \f]
0080         because the difference
0081         \f[
0082             \left(\frac{1}{n-1} - \frac{1}{n}\right)(x_n - \mu_n)^2 = \frac{1}{n(n-1)}(x_n - \mu_n)^2.
0083         \f]
0084         converges to zero as \f$n \rightarrow \infty\f$. However, for small \f$ n \f$ the difference
0085         can be non-negligible.
0086     */
0087     template<typename Sample, typename MeanFeature, typename Tag>
0088     struct variance_impl
0089       : accumulator_base
0090     {
0091         // for boost::result_of
0092         typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type result_type;
0093
0094         template<typename Args>
0095         variance_impl(Args const &args)
0096           : variance(numeric::fdiv(args[sample | Sample()], numeric::one<std::size_t>::value))
0097         {
0098         }
0099
0100         template<typename Args>
0101         void operator ()(Args const &args)
0102         {
0103             std::size_t cnt = count(args);
0104
0105             if(cnt > 1)
0106             {
0107                 extractor<MeanFeature> mean;
0108                 result_type tmp = args[parameter::keyword<Tag>::get()] - mean(args);
0109                 this->variance =
0110                     numeric::fdiv(this->variance * (cnt - 1), cnt)
0111                   + numeric::fdiv(tmp * tmp, cnt - 1);
0112             }
0113         }
0114
0115         result_type result(dont_care) const
0116         {
0117             return this->variance;
0118         }
0119
0120         // make this accumulator serializeable
0121         template<class Archive>
0122         void serialize(Archive & ar, const unsigned int /* file_version */)
0123         {
0124             ar & variance;
0125         }
0126
0127     private:
0128         result_type variance;
0129     };
0130
0131 } // namespace impl
0132
0133 ///////////////////////////////////////////////////////////////////////////////
0134 // tag::variance
0135 // tag::immediate_variance
0136 //
0137 namespace tag
0138 {
0139     struct lazy_variance
0140       : depends_on<moment<2>, mean>
0141     {
0142         /// INTERNAL ONLY
0143         ///
0144         typedef accumulators::impl::lazy_variance_impl<mpl::_1, mean> impl;
0145     };
0146
0147     struct variance
0148       : depends_on<count, immediate_mean>
0149     {
0150         /// INTERNAL ONLY
0151         ///
0152         typedef accumulators::impl::variance_impl<mpl::_1, mean, sample> impl;
0153     };
0154 }
0155
0156 ///////////////////////////////////////////////////////////////////////////////
0157 // extract::lazy_variance
0158 // extract::variance
0159 //
0160 namespace extract
0161 {
0162     extractor<tag::lazy_variance> const lazy_variance = {};
0163     extractor<tag::variance> const variance = {};
0164
0165     BOOST_ACCUMULATORS_IGNORE_GLOBAL(lazy_variance)
0166     BOOST_ACCUMULATORS_IGNORE_GLOBAL(variance)
0167 }
0168
0169 using extract::lazy_variance;
0170 using extract::variance;
0171
0172 // variance(lazy) -> lazy_variance
0173 template<>
0174 struct as_feature<tag::variance(lazy)>
0175 {
0176     typedef tag::lazy_variance type;
0177 };
0178
0179 // variance(immediate) -> variance
0180 template<>
0181 struct as_feature<tag::variance(immediate)>
0182 {
0183     typedef tag::variance type;
0184 };
0185
0186 // for the purposes of feature-based dependency resolution,
0187 // immediate_variance provides the same feature as variance
0188 template<>
0189 struct feature_of<tag::lazy_variance>
0190   : feature_of<tag::variance>
0191 {
0192 };
0193
0194 // So that variance can be automatically substituted with
0195 // weighted_variance when the weight parameter is non-void.
0196 template<>
0197 struct as_weighted_feature<tag::variance>
0198 {
0199     typedef tag::weighted_variance type;
0200 };
0201
0202 // for the purposes of feature-based dependency resolution,
0203 // weighted_variance provides the same feature as variance
0204 template<>
0205 struct feature_of<tag::weighted_variance>
0206   : feature_of<tag::variance>
0207 {
0208 };
0209
0210 // So that immediate_variance can be automatically substituted with
0211 // immediate_weighted_variance when the weight parameter is non-void.
0212 template<>
0213 struct as_weighted_feature<tag::lazy_variance>
0214 {
0215     typedef tag::lazy_weighted_variance type;
0216 };
0217
0218 // for the purposes of feature-based dependency resolution,
0219 // immediate_weighted_variance provides the same feature as immediate_variance
0220 template<>
0221 struct feature_of<tag::lazy_weighted_variance>
0222   : feature_of<tag::lazy_variance>
0223 {
0224 };
0225
0226 ////////////////////////////////////////////////////////////////////////////
0227 //// droppable_accumulator<variance_impl>
0228 ////  need to specialize droppable lazy variance to cache the result at the
0229 ////  point the accumulator is dropped.
0230 ///// INTERNAL ONLY
0231 /////
0232 //template<typename Sample, typename MeanFeature>
0233 //struct droppable_accumulator<impl::variance_impl<Sample, MeanFeature> >
0234 //  : droppable_accumulator_base<
0235 //        with_cached_result<impl::variance_impl<Sample, MeanFeature> >
0236 //    >
0237 //{
0238 //    template<typename Args>
0239 //    droppable_accumulator(Args const &args)
0240 //      : droppable_accumulator::base(args)
0241 //    {
0242 //    }
0243 //};
0244
0245 }} // namespace boost::accumulators
0246
0247 #endif