1 /////////////////////////////////////////////////////////////////////////////// 2 // rolling_variance.hpp 3 // Copyright (C) 2005 Eric Niebler 4 // Copyright (C) 2014 Pieter Bastiaan Ober (Integricom). 5 // Distributed under the Boost Software License, Version 1.0. 6 // (See accompanying file LICENSE_1_0.txt or copy at 7 // http://www.boost.org/LICENSE_1_0.txt) 8 9 #ifndef BOOST_ACCUMULATORS_STATISTICS_ROLLING_VARIANCE_HPP_EAN_15_11_2011 10 #define BOOST_ACCUMULATORS_STATISTICS_ROLLING_VARIANCE_HPP_EAN_15_11_2011 11 12 #include <boost/accumulators/accumulators.hpp> 13 #include <boost/accumulators/statistics/stats.hpp> 14 15 #include <boost/mpl/placeholders.hpp> 16 #include <boost/accumulators/framework/accumulator_base.hpp> 17 #include <boost/accumulators/framework/extractor.hpp> 18 #include <boost/accumulators/numeric/functional.hpp> 19 #include <boost/accumulators/framework/parameters/sample.hpp> 20 #include <boost/accumulators/framework/depends_on.hpp> 21 #include <boost/accumulators/statistics_fwd.hpp> 22 #include <boost/accumulators/statistics/rolling_mean.hpp> 23 #include <boost/accumulators/statistics/rolling_moment.hpp> 24 25 #include <boost/type_traits/is_arithmetic.hpp> 26 #include <boost/utility/enable_if.hpp> 27 28 namespace boost { namespace accumulators 29 { 30 namespace impl 31 { 32 //! Immediate (lazy) calculation of the rolling variance. 33 /*! 34 Calculation of sample variance \f$\sigma_n^2\f$ is done as follows, see also 35 http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance. 36 For a rolling window of size \f$N\f$, when \f$n <= N\f$, the variance is computed according to the formula 37 \f[ 38 \sigma_n^2 = \frac{1}{n-1} \sum_{i = 1}^n (x_i - \mu_n)^2. 39 \f] 40 When \f$n > N\f$, the sample variance over the window becomes: 41 \f[ 42 \sigma_n^2 = \frac{1}{N-1} \sum_{i = n-N+1}^n (x_i - \mu_n)^2. 43 \f] 44 */ 45 /////////////////////////////////////////////////////////////////////////////// 46 // lazy_rolling_variance_impl 47 // 48 template<typename Sample> 49 struct lazy_rolling_variance_impl 50 : accumulator_base 51 { 52 // for boost::result_of 53 typedef typename numeric::functional::fdiv<Sample, std::size_t,void,void>::result_type result_type; 54 lazy_rolling_variance_implboost::accumulators::impl::lazy_rolling_variance_impl55 lazy_rolling_variance_impl(dont_care) {} 56 57 template<typename Args> resultboost::accumulators::impl::lazy_rolling_variance_impl58 result_type result(Args const &args) const 59 { 60 result_type mean = rolling_mean(args); 61 size_t nr_samples = rolling_count(args); 62 if (nr_samples < 2) return result_type(); 63 return nr_samples*(rolling_moment<2>(args) - mean*mean)/(nr_samples-1); 64 } 65 66 // serialization is done by accumulators it depends on 67 template<class Archive> serializeboost::accumulators::impl::lazy_rolling_variance_impl68 void serialize(Archive & ar, const unsigned int file_version) {} 69 }; 70 71 //! Iterative calculation of the rolling variance. 72 /*! 73 Iterative calculation of sample variance \f$\sigma_n^2\f$ is done as follows, see also 74 http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance. 75 For a rolling window of size \f$N\f$, for the first \f$N\f$ samples, the variance is computed according to the formula 76 \f[ 77 \sigma_n^2 = \frac{1}{n-1} \sum_{i = 1}^n (x_i - \mu_n)^2 = \frac{1}{n-1}M_{2,n}, 78 \f] 79 where the sum of squares \f$M_{2,n}\f$ can be recursively computed as: 80 \f[ 81 M_{2,n} = \sum_{i = 1}^n (x_i - \mu_n)^2 = M_{2,n-1} + (x_n - \mu_n)(x_n - \mu_{n-1}), 82 \f] 83 and the estimate of the sample mean as: 84 \f[ 85 \mu_n = \frac{1}{n} \sum_{i = 1}^n x_i = \mu_{n-1} + \frac{1}{n}(x_n - \mu_{n-1}). 86 \f] 87 For further samples, when the rolling window is fully filled with data, one has to take into account that the oldest 88 sample \f$x_{n-N}\f$ is dropped from the window. The sample variance over the window now becomes: 89 \f[ 90 \sigma_n^2 = \frac{1}{N-1} \sum_{i = n-N+1}^n (x_i - \mu_n)^2 = \frac{1}{n-1}M_{2,n}, 91 \f] 92 where the sum of squares \f$M_{2,n}\f$ now equals: 93 \f[ 94 M_{2,n} = \sum_{i = n-N+1}^n (x_i - \mu_n)^2 = M_{2,n-1} + (x_n - \mu_n)(x_n - \mu_{n-1}) - (x_{n-N} - \mu_n)(x_{n-N} - \mu_{n-1}), 95 \f] 96 and the estimated mean is: 97 \f[ 98 \mu_n = \frac{1}{N} \sum_{i = n-N+1}^n x_i = \mu_{n-1} + \frac{1}{n}(x_n - x_{n-N}). 99 \f] 100 101 Note that the sample variance is not defined for \f$n <= 1\f$. 102 103 */ 104 /////////////////////////////////////////////////////////////////////////////// 105 // immediate_rolling_variance_impl 106 // 107 template<typename Sample> 108 struct immediate_rolling_variance_impl 109 : accumulator_base 110 { 111 // for boost::result_of 112 typedef typename numeric::functional::fdiv<Sample, std::size_t>::result_type result_type; 113 114 template<typename Args> immediate_rolling_variance_implboost::accumulators::impl::immediate_rolling_variance_impl115 immediate_rolling_variance_impl(Args const &args) 116 : previous_mean_(numeric::fdiv(args[sample | Sample()], numeric::one<std::size_t>::value)) 117 , sum_of_squares_(numeric::fdiv(args[sample | Sample()], numeric::one<std::size_t>::value)) 118 { 119 } 120 121 template<typename Args> operator ()boost::accumulators::impl::immediate_rolling_variance_impl122 void operator()(Args const &args) 123 { 124 Sample added_sample = args[sample]; 125 126 result_type mean = immediate_rolling_mean(args); 127 sum_of_squares_ += (added_sample-mean)*(added_sample-previous_mean_); 128 129 if(is_rolling_window_plus1_full(args)) 130 { 131 Sample removed_sample = rolling_window_plus1(args).front(); 132 sum_of_squares_ -= (removed_sample-mean)*(removed_sample-previous_mean_); 133 prevent_underflow(sum_of_squares_); 134 } 135 previous_mean_ = mean; 136 } 137 138 template<typename Args> resultboost::accumulators::impl::immediate_rolling_variance_impl139 result_type result(Args const &args) const 140 { 141 size_t nr_samples = rolling_count(args); 142 if (nr_samples < 2) return result_type(); 143 return numeric::fdiv(sum_of_squares_,(nr_samples-1)); 144 } 145 146 // make this accumulator serializeable 147 template<class Archive> serializeboost::accumulators::impl::immediate_rolling_variance_impl148 void serialize(Archive & ar, const unsigned int file_version) 149 { 150 ar & previous_mean_; 151 ar & sum_of_squares_; 152 } 153 154 private: 155 156 result_type previous_mean_; 157 result_type sum_of_squares_; 158 159 template<typename T> prevent_underflowboost::accumulators::impl::immediate_rolling_variance_impl160 void prevent_underflow(T &non_negative_number,typename boost::enable_if<boost::is_arithmetic<T>,T>::type* = 0) 161 { 162 if (non_negative_number < T(0)) non_negative_number = T(0); 163 } 164 template<typename T> prevent_underflowboost::accumulators::impl::immediate_rolling_variance_impl165 void prevent_underflow(T &non_arithmetic_quantity,typename boost::disable_if<boost::is_arithmetic<T>,T>::type* = 0) 166 { 167 } 168 }; 169 } // namespace impl 170 171 /////////////////////////////////////////////////////////////////////////////// 172 // tag:: lazy_rolling_variance 173 // tag:: immediate_rolling_variance 174 // tag:: rolling_variance 175 // 176 namespace tag 177 { 178 struct lazy_rolling_variance 179 : depends_on< rolling_count, rolling_mean, rolling_moment<2> > 180 { 181 /// INTERNAL ONLY 182 /// 183 typedef accumulators::impl::lazy_rolling_variance_impl< mpl::_1 > impl; 184 185 #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED 186 /// tag::rolling_window::window_size named parameter 187 static boost::parameter::keyword<tag::rolling_window_size> const window_size; 188 #endif 189 }; 190 191 struct immediate_rolling_variance 192 : depends_on< rolling_window_plus1, rolling_count, immediate_rolling_mean> 193 { 194 /// INTERNAL ONLY 195 /// 196 typedef accumulators::impl::immediate_rolling_variance_impl< mpl::_1> impl; 197 198 #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED 199 /// tag::rolling_window::window_size named parameter 200 static boost::parameter::keyword<tag::rolling_window_size> const window_size; 201 #endif 202 }; 203 204 // make immediate_rolling_variance the default implementation 205 struct rolling_variance : immediate_rolling_variance {}; 206 } // namespace tag 207 208 /////////////////////////////////////////////////////////////////////////////// 209 // extract::lazy_rolling_variance 210 // extract::immediate_rolling_variance 211 // extract::rolling_variance 212 // 213 namespace extract 214 { 215 extractor<tag::lazy_rolling_variance> const lazy_rolling_variance = {}; 216 extractor<tag::immediate_rolling_variance> const immediate_rolling_variance = {}; 217 extractor<tag::rolling_variance> const rolling_variance = {}; 218 219 BOOST_ACCUMULATORS_IGNORE_GLOBAL(lazy_rolling_variance) 220 BOOST_ACCUMULATORS_IGNORE_GLOBAL(immediate_rolling_variance) 221 BOOST_ACCUMULATORS_IGNORE_GLOBAL(rolling_variance) 222 } 223 224 using extract::lazy_rolling_variance; 225 using extract::immediate_rolling_variance; 226 using extract::rolling_variance; 227 228 // rolling_variance(lazy) -> lazy_rolling_variance 229 template<> 230 struct as_feature<tag::rolling_variance(lazy)> 231 { 232 typedef tag::lazy_rolling_variance type; 233 }; 234 235 // rolling_variance(immediate) -> immediate_rolling_variance 236 template<> 237 struct as_feature<tag::rolling_variance(immediate)> 238 { 239 typedef tag::immediate_rolling_variance type; 240 }; 241 242 // for the purposes of feature-based dependency resolution, 243 // lazy_rolling_variance provides the same feature as rolling_variance 244 template<> 245 struct feature_of<tag::lazy_rolling_variance> 246 : feature_of<tag::rolling_variance> 247 { 248 }; 249 250 // for the purposes of feature-based dependency resolution, 251 // immediate_rolling_variance provides the same feature as rolling_variance 252 template<> 253 struct feature_of<tag::immediate_rolling_variance> 254 : feature_of<tag::rolling_variance> 255 { 256 }; 257 }} // namespace boost::accumulators 258 259 #endif 260