]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | : std::binary_function< | |
54 | Left | |
55 | , Right | |
56 | , ublas::matrix< | |
57 | typename functional::multiplies< | |
58 | typename Left::value_type | |
59 | , typename Right::value_type | |
60 | >::result_type | |
61 | > | |
62 | > | |
63 | { | |
64 | typedef | |
65 | ublas::matrix< | |
66 | typename functional::multiplies< | |
67 | typename Left::value_type | |
68 | , typename Right::value_type | |
69 | >::result_type | |
70 | > | |
71 | result_type; | |
72 | ||
73 | result_type | |
74 | operator ()(Left & left, Right & right) const | |
75 | { | |
76 | std::size_t left_size = left.size(); | |
77 | std::size_t right_size = right.size(); | |
78 | result_type result(left_size, right_size); | |
79 | for (std::size_t i = 0; i < left_size; ++i) | |
80 | for (std::size_t j = 0; j < right_size; ++j) | |
81 | result(i,j) = numeric::multiplies(left[i], right[j]); | |
82 | return result; | |
83 | } | |
84 | }; | |
85 | } | |
86 | ||
87 | namespace op | |
88 | { | |
89 | struct outer_product | |
90 | : boost::detail::function2<functional::outer_product<_1, _2, functional::tag<_1>, functional::tag<_2> > > | |
91 | {}; | |
92 | } | |
93 | ||
94 | namespace | |
95 | { | |
96 | op::outer_product const &outer_product = boost::detail::pod_singleton<op::outer_product>::instance; | |
97 | } | |
98 | ||
99 | }} | |
100 | ||
101 | namespace boost { namespace accumulators | |
102 | { | |
103 | ||
104 | namespace impl | |
105 | { | |
106 | /////////////////////////////////////////////////////////////////////////////// | |
107 | // covariance_impl | |
108 | // | |
109 | /** | |
110 | @brief Covariance Estimator | |
111 | ||
112 | An iterative Monte Carlo estimator for the covariance \f$\mathrm{Cov}(X,X')\f$, where \f$X\f$ is a sample | |
113 | and \f$X'\f$ is a variate, is given by: | |
114 | ||
115 | \f[ | |
116 | \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, | |
117 | \f] | |
118 | ||
119 | \f$\hat{\mu}_n\f$ and \f$\hat{\mu}_n'\f$ being the means of the samples and variates. | |
120 | */ | |
121 | template<typename Sample, typename VariateType, typename VariateTag> | |
122 | struct covariance_impl | |
123 | : accumulator_base | |
124 | { | |
125 | typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type sample_type; | |
126 | typedef typename numeric::functional::fdiv<VariateType, std::size_t>::result_type variate_type; | |
127 | // for boost::result_of | |
128 | typedef typename numeric::functional::outer_product<sample_type, variate_type>::result_type result_type; | |
129 | ||
130 | template<typename Args> | |
131 | covariance_impl(Args const &args) | |
132 | : cov_( | |
133 | numeric::outer_product( | |
134 | numeric::fdiv(args[sample | Sample()], (std::size_t)1) | |
135 | , numeric::fdiv(args[parameter::keyword<VariateTag>::get() | VariateType()], (std::size_t)1) | |
136 | ) | |
137 | ) | |
138 | { | |
139 | } | |
140 | ||
141 | template<typename Args> | |
142 | void operator ()(Args const &args) | |
143 | { | |
144 | std::size_t cnt = count(args); | |
145 | ||
146 | if (cnt > 1) | |
147 | { | |
148 | extractor<tag::mean_of_variates<VariateType, VariateTag> > const some_mean_of_variates = {}; | |
149 | ||
150 | this->cov_ = this->cov_*(cnt-1.)/cnt | |
151 | + numeric::outer_product( | |
152 | some_mean_of_variates(args) - args[parameter::keyword<VariateTag>::get()] | |
153 | , mean(args) - args[sample] | |
154 | ) / (cnt-1.); | |
155 | } | |
156 | } | |
157 | ||
158 | result_type result(dont_care) const | |
159 | { | |
160 | return this->cov_; | |
161 | } | |
162 | ||
163 | private: | |
164 | result_type cov_; | |
165 | }; | |
166 | ||
167 | } // namespace impl | |
168 | ||
169 | /////////////////////////////////////////////////////////////////////////////// | |
170 | // tag::covariance | |
171 | // | |
172 | namespace tag | |
173 | { | |
174 | template<typename VariateType, typename VariateTag> | |
175 | struct covariance | |
176 | : depends_on<count, mean, mean_of_variates<VariateType, VariateTag> > | |
177 | { | |
178 | typedef accumulators::impl::covariance_impl<mpl::_1, VariateType, VariateTag> impl; | |
179 | }; | |
180 | ||
181 | struct abstract_covariance | |
182 | : depends_on<> | |
183 | { | |
184 | }; | |
185 | } | |
186 | ||
187 | /////////////////////////////////////////////////////////////////////////////// | |
188 | // extract::covariance | |
189 | // | |
190 | namespace extract | |
191 | { | |
192 | extractor<tag::abstract_covariance> const covariance = {}; | |
193 | ||
194 | BOOST_ACCUMULATORS_IGNORE_GLOBAL(covariance) | |
195 | } | |
196 | ||
197 | using extract::covariance; | |
198 | ||
199 | template<typename VariateType, typename VariateTag> | |
200 | struct feature_of<tag::covariance<VariateType, VariateTag> > | |
201 | : feature_of<tag::abstract_covariance> | |
202 | { | |
203 | }; | |
204 | ||
205 | // So that covariance can be automatically substituted with | |
206 | // weighted_covariance when the weight parameter is non-void. | |
207 | template<typename VariateType, typename VariateTag> | |
208 | struct as_weighted_feature<tag::covariance<VariateType, VariateTag> > | |
209 | { | |
210 | typedef tag::weighted_covariance<VariateType, VariateTag> type; | |
211 | }; | |
212 | ||
213 | template<typename VariateType, typename VariateTag> | |
214 | struct feature_of<tag::weighted_covariance<VariateType, VariateTag> > | |
215 | : feature_of<tag::covariance<VariateType, VariateTag> > | |
216 | {}; | |
217 | ||
218 | }} // namespace boost::accumulators | |
219 | ||
220 | #endif |