accumulators/statistics/weighted_covariance.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // weighted_covariance.hpp
0003 //
0004 //  Copyright 2006 Daniel Egloff, Olivier Gygi. 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_WEIGHTED_COVARIANCE_HPP_DE_01_01_2006
0009 #define BOOST_ACCUMULATORS_STATISTICS_WEIGHTED_COVARIANCE_HPP_DE_01_01_2006
0010
0011 #include <vector>
0012 #include <limits>
0013 #include <numeric>
0014 #include <functional>
0015 #include <complex>
0016 #include <boost/mpl/assert.hpp>
0017 #include <boost/mpl/bool.hpp>
0018 #include <boost/range.hpp>
0019 #include <boost/parameter/keyword.hpp>
0020 #include <boost/mpl/placeholders.hpp>
0021 #include <boost/numeric/ublas/io.hpp>
0022 #include <boost/numeric/ublas/matrix.hpp>
0023 #include <boost/type_traits/is_scalar.hpp>
0024 #include <boost/type_traits/is_same.hpp>
0025 #include <boost/accumulators/framework/accumulator_base.hpp>
0026 #include <boost/accumulators/framework/extractor.hpp>
0027 #include <boost/accumulators/numeric/functional.hpp>
0028 #include <boost/accumulators/framework/parameters/sample.hpp>
0029 #include <boost/accumulators/statistics_fwd.hpp>
0030 #include <boost/accumulators/statistics/count.hpp>
0031 #include <boost/accumulators/statistics/covariance.hpp> // for numeric::outer_product() and type traits
0032 #include <boost/accumulators/statistics/weighted_mean.hpp>
0033
0034 namespace boost { namespace accumulators
0035 {
0036
0037 namespace impl
0038 {
0039     ///////////////////////////////////////////////////////////////////////////////
0040     // weighted_covariance_impl
0041     //
0042     /**
0043         @brief Weighted Covariance Estimator
0044
0045         An iterative Monte Carlo estimator for the weighted covariance \f$\mathrm{Cov}(X,X')\f$, where \f$X\f$ is a sample
0046         and \f$X'\f$ a variate, is given by:
0047
0048         \f[
0049             \hat{c}_n = \frac{\bar{w}_n-w_n}{\bar{w}_n} \hat{c}_{n-1} + \frac{w_n}{\bar{w}_n-w_n}(X_n - \hat{\mu}_n)(X_n' - \hat{\mu}_n'),
0050             \quad n\ge2,\quad\hat{c}_1 = 0,
0051         \f]
0052
0053         \f$\hat{\mu}_n\f$ and \f$\hat{\mu}_n'\f$ being the weighted means of the samples and variates and
0054         \f$\bar{w}_n\f$ the sum of the \f$n\f$ first weights \f$w_i\f$.
0055     */
0056     template<typename Sample, typename Weight, typename VariateType, typename VariateTag>
0057     struct weighted_covariance_impl
0058       : accumulator_base
0059     {
0060         typedef typename numeric::functional::multiplies<Weight, typename numeric::functional::fdiv<Sample, std::size_t>::result_type>::result_type weighted_sample_type;
0061         typedef typename numeric::functional::multiplies<Weight, typename numeric::functional::fdiv<VariateType, std::size_t>::result_type>::result_type weighted_variate_type;
0062         // for boost::result_of
0063         typedef typename numeric::functional::outer_product<weighted_sample_type, weighted_variate_type>::result_type result_type;
0064
0065         template<typename Args>
0066         weighted_covariance_impl(Args const &args)
0067           : cov_(
0068                 numeric::outer_product(
0069                     numeric::fdiv(args[sample | Sample()], (std::size_t)1)
0070                       * numeric::one<Weight>::value
0071                   , numeric::fdiv(args[parameter::keyword<VariateTag>::get() | VariateType()], (std::size_t)1)
0072                       * numeric::one<Weight>::value
0073                 )
0074             )
0075         {
0076         }
0077
0078         template<typename Args>
0079         void operator ()(Args const &args)
0080         {
0081             std::size_t cnt = count(args);
0082
0083             if (cnt > 1)
0084             {
0085                 extractor<tag::weighted_mean_of_variates<VariateType, VariateTag> > const some_weighted_mean_of_variates = {};
0086
0087                 this->cov_ = this->cov_ * (sum_of_weights(args) - args[weight]) / sum_of_weights(args)
0088                            + numeric::outer_product(
0089                                  some_weighted_mean_of_variates(args) - args[parameter::keyword<VariateTag>::get()]
0090                                , weighted_mean(args) - args[sample]
0091                              ) * args[weight] / (sum_of_weights(args) - args[weight]);
0092             }
0093         }
0094
0095         result_type result(dont_care) const
0096         {
0097             return this->cov_;
0098         }
0099
0100         // make this accumulator serializeable
0101         template<class Archive>
0102         void serialize(Archive & ar, const unsigned int file_version)
0103         {
0104             ar & cov_;
0105         }
0106
0107     private:
0108         result_type cov_;
0109     };
0110
0111 } // namespace impl
0112
0113 ///////////////////////////////////////////////////////////////////////////////
0114 // tag::weighted_covariance
0115 //
0116 namespace tag
0117 {
0118     template<typename VariateType, typename VariateTag>
0119     struct weighted_covariance
0120       : depends_on<count, sum_of_weights, weighted_mean, weighted_mean_of_variates<VariateType, VariateTag> >
0121     {
0122         typedef accumulators::impl::weighted_covariance_impl<mpl::_1, mpl::_2, VariateType, VariateTag> impl;
0123     };
0124 }
0125
0126 ///////////////////////////////////////////////////////////////////////////////
0127 // extract::weighted_covariance
0128 //
0129 namespace extract
0130 {
0131     extractor<tag::abstract_covariance> const weighted_covariance = {};
0132
0133     BOOST_ACCUMULATORS_IGNORE_GLOBAL(weighted_covariance)
0134 }
0135
0136 using extract::weighted_covariance;
0137
0138 }} // namespace boost::accumulators
0139
0140 #endif