1 /////////////////////////////////////////////////////////////////////////////// 2 // covariance.hpp 3 // 4 // Copyright 2006 Daniel Egloff, Olivier Gygi. 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_COVARIANCE_HPP_DE_01_01_2006 9 #define BOOST_ACCUMULATORS_STATISTICS_COVARIANCE_HPP_DE_01_01_2006 10 11 #include <vector> 12 #include <limits> 13 #include <numeric> 14 #include <functional> 15 #include <complex> 16 #include <boost/mpl/assert.hpp> 17 #include <boost/mpl/bool.hpp> 18 #include <boost/range.hpp> 19 #include <boost/parameter/keyword.hpp> 20 #include <boost/mpl/placeholders.hpp> 21 #include <boost/numeric/ublas/io.hpp> 22 #include <boost/numeric/ublas/matrix.hpp> 23 #include <boost/type_traits/is_scalar.hpp> 24 #include <boost/type_traits/is_same.hpp> 25 #include <boost/accumulators/framework/accumulator_base.hpp> 26 #include <boost/accumulators/framework/extractor.hpp> 27 #include <boost/accumulators/numeric/functional.hpp> 28 #include <boost/accumulators/framework/parameters/sample.hpp> 29 #include <boost/accumulators/statistics_fwd.hpp> 30 #include <boost/accumulators/statistics/count.hpp> 31 #include <boost/accumulators/statistics/mean.hpp> 32 33 namespace boost { namespace numeric 34 { 35 namespace functional 36 { 37 struct std_vector_tag; 38 39 /////////////////////////////////////////////////////////////////////////////// 40 // functional::outer_product 41 template<typename Left, typename Right, typename EnableIf = void> 42 struct outer_product_base 43 : functional::multiplies<Left, Right> 44 {}; 45 46 template<typename Left, typename Right, typename LeftTag = typename tag<Left>::type, typename RightTag = typename tag<Right>::type> 47 struct outer_product 48 : outer_product_base<Left, Right, void> 49 {}; 50 51 template<typename Left, typename Right> 52 struct outer_product<Left, Right, std_vector_tag, std_vector_tag> 53 { 54 typedef Left first_argument_type; 55 typedef Right second_argument_type; 56 typedef 57 ublas::matrix< 58 typename functional::multiplies< 59 typename Left::value_type 60 , typename Right::value_type 61 >::result_type 62 > 63 result_type; 64 65 result_type operator ()boost::numeric::functional::outer_product66 operator ()(Left & left, Right & right) const 67 { 68 std::size_t left_size = left.size(); 69 std::size_t right_size = right.size(); 70 result_type result(left_size, right_size); 71 for (std::size_t i = 0; i < left_size; ++i) 72 for (std::size_t j = 0; j < right_size; ++j) 73 result(i,j) = numeric::multiplies(left[i], right[j]); 74 return result; 75 } 76 }; 77 } 78 79 namespace op 80 { 81 struct outer_product 82 : boost::detail::function2<functional::outer_product<_1, _2, functional::tag<_1>, functional::tag<_2> > > 83 {}; 84 } 85 86 namespace 87 { 88 op::outer_product const &outer_product = boost::detail::pod_singleton<op::outer_product>::instance; 89 } 90 91 }} 92 93 namespace boost { namespace accumulators 94 { 95 96 namespace impl 97 { 98 /////////////////////////////////////////////////////////////////////////////// 99 // covariance_impl 100 // 101 /** 102 @brief Covariance Estimator 103 104 An iterative Monte Carlo estimator for the covariance \f$\mathrm{Cov}(X,X')\f$, where \f$X\f$ is a sample 105 and \f$X'\f$ is a variate, is given by: 106 107 \f[ 108 \hat{c}_n = \frac{n-1}{n} \hat{c}_{n-1} + \frac{1}{n-1}(X_n - \hat{\mu}_n)(X_n' - \hat{\mu}_n'),\quad n\ge2,\quad\hat{c}_1 = 0, 109 \f] 110 111 \f$\hat{\mu}_n\f$ and \f$\hat{\mu}_n'\f$ being the means of the samples and variates. 112 */ 113 template<typename Sample, typename VariateType, typename VariateTag> 114 struct covariance_impl 115 : accumulator_base 116 { 117 typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type sample_type; 118 typedef typename numeric::functional::fdiv<VariateType, std::size_t>::result_type variate_type; 119 // for boost::result_of 120 typedef typename numeric::functional::outer_product<sample_type, variate_type>::result_type result_type; 121 122 template<typename Args> covariance_implboost::accumulators::impl::covariance_impl123 covariance_impl(Args const &args) 124 : cov_( 125 numeric::outer_product( 126 numeric::fdiv(args[sample | Sample()], (std::size_t)1) 127 , numeric::fdiv(args[parameter::keyword<VariateTag>::get() | VariateType()], (std::size_t)1) 128 ) 129 ) 130 { 131 } 132 133 template<typename Args> operator ()boost::accumulators::impl::covariance_impl134 void operator ()(Args const &args) 135 { 136 std::size_t cnt = count(args); 137 138 if (cnt > 1) 139 { 140 extractor<tag::mean_of_variates<VariateType, VariateTag> > const some_mean_of_variates = {}; 141 142 this->cov_ = this->cov_*(cnt-1.)/cnt 143 + numeric::outer_product( 144 some_mean_of_variates(args) - args[parameter::keyword<VariateTag>::get()] 145 , mean(args) - args[sample] 146 ) / (cnt-1.); 147 } 148 } 149 resultboost::accumulators::impl::covariance_impl150 result_type result(dont_care) const 151 { 152 return this->cov_; 153 } 154 155 // make this accumulator serializeable 156 template<class Archive> serializeboost::accumulators::impl::covariance_impl157 void serialize(Archive & ar, const unsigned int file_version) 158 { 159 ar & cov_; 160 } 161 162 private: 163 result_type cov_; 164 }; 165 166 } // namespace impl 167 168 /////////////////////////////////////////////////////////////////////////////// 169 // tag::covariance 170 // 171 namespace tag 172 { 173 template<typename VariateType, typename VariateTag> 174 struct covariance 175 : depends_on<count, mean, mean_of_variates<VariateType, VariateTag> > 176 { 177 typedef accumulators::impl::covariance_impl<mpl::_1, VariateType, VariateTag> impl; 178 }; 179 180 struct abstract_covariance 181 : depends_on<> 182 { 183 }; 184 } 185 186 /////////////////////////////////////////////////////////////////////////////// 187 // extract::covariance 188 // 189 namespace extract 190 { 191 extractor<tag::abstract_covariance> const covariance = {}; 192 193 BOOST_ACCUMULATORS_IGNORE_GLOBAL(covariance) 194 } 195 196 using extract::covariance; 197 198 template<typename VariateType, typename VariateTag> 199 struct feature_of<tag::covariance<VariateType, VariateTag> > 200 : feature_of<tag::abstract_covariance> 201 { 202 }; 203 204 // So that covariance can be automatically substituted with 205 // weighted_covariance when the weight parameter is non-void. 206 template<typename VariateType, typename VariateTag> 207 struct as_weighted_feature<tag::covariance<VariateType, VariateTag> > 208 { 209 typedef tag::weighted_covariance<VariateType, VariateTag> type; 210 }; 211 212 template<typename VariateType, typename VariateTag> 213 struct feature_of<tag::weighted_covariance<VariateType, VariateTag> > 214 : feature_of<tag::covariance<VariateType, VariateTag> > 215 {}; 216 217 }} // namespace boost::accumulators 218 219 #endif 220