Ticket #6128: RollingVariance.hpp

///////////////////////////////////////////////////////////////////////////////
// rolling_variance.hpp
//  Copyright (C) 2005 Eric Niebler    
//  Copyright (C) 2011 Pieter Bastiaan Ober (Integricom). Distributed under the Boost
//  Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_ACCUMULATORS_STATISTICS_ROLLING_VARIANCE_HPP_EAN_15_11_2011
#define BOOST_ACCUMULATORS_STATISTICS_ROLLING_VARIANCE_HPP_EAN_15_11_2011

#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>

#include <boost/mpl/placeholders.hpp>
#include <boost/accumulators/framework/accumulator_base.hpp>
#include <boost/accumulators/framework/extractor.hpp>
#include <boost/accumulators/numeric/functional.hpp>
#include <boost/accumulators/framework/parameters/sample.hpp>
#include <boost/accumulators/framework/depends_on.hpp>
#include <boost/accumulators/statistics_fwd.hpp>
#include <boost/accumulators/statistics/rolling_window.hpp>
#include <boost/accumulators/statistics/rolling_count.hpp>

#include <RollingMean.hpp>

// See: Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1983). 
// Algorithms for Computing the Sample Variance: Analysis and Recommendations. 
// The American Statistician 37, 242-247.
// 
// Variance(x[1:N]):
// {
//    M[1] = x[1]
//    S[1] = 0.0
// 
//              j=2:N
//              {
//                      M[j] = M[j-1] + (1/j)*(x[j]-M[j]-1)
//                      S[j] = S[j-1] + ((j-1)/j)*(x[j]-M[j-1])^2
//                      Variance = S[j]/j or
//                      Variance = S[j]/(j-1)
//              }
// }
//    
// When the buffer of size N is full, not only points are added but points are dropped as well.
// from the window over which the variance is computed. Adaptation to a rolling window gives
// for j=N+1,...: 
//      { 
//                      M[j] = M[j-1] + (1/N)*(x[j]-x[j-N])
//                      S[j] = S[j-1] + (N/(N+1))*(x[j]-M[j-1])^2 - (N/(N+1))*(x[j-N]-M[j])^2
//                      Variance = S[j]/N or
//                      Variance = S[j]/(N-1)
// }


namespace boost { namespace accumulators
{
namespace impl
{
    ///////////////////////////////////////////////////////////////////////////////
    // rolling_variance_impl
    // 
    template<typename Sample>
    struct rolling_variance_impl
      : accumulator_base
    {
        typedef Sample result_type;

        template<typename Args>
        rolling_variance_impl(Args const &args)
          : previous_mean_(0.0)
                         ,      sum_of_squares_(0.0)
        {
                          //    VOID;
                  }

        template<typename Args>
        void operator()(Args const &args)
        {
                                Sample added_sample = args[sample];
                                size_t nr_samples = rolling_count(args);

                                Sample mean = rolling_mean(args);

            if(is_rolling_window_plus1_full(args))
            {
                                        Sample weight = static_cast<Sample>(nr_samples)/static_cast<Sample>(nr_samples+1);
                                        Sample removed_sample = rolling_window_plus1(args).front();
                                        sum_of_squares_ += weight*(
                                                                                                (added_sample-previous_mean_)*(added_sample-previous_mean_) -
                                                                                                (removed_sample-mean)*(removed_sample-mean)
                                                                                                );
                                        if (sum_of_squares_ < 0.0)      sum_of_squares_ = 0.0;
            }
                                else
                                {
                                        Sample weight = static_cast<Sample>(nr_samples-1)/static_cast<Sample>(nr_samples);
                                        sum_of_squares_ += weight*(added_sample-previous_mean_)*(added_sample-previous_mean_);
                                }
                                previous_mean_ = mean;
        }

        template<typename Args>
        result_type result(Args const &args) const
        {
                          size_t nr_samples = rolling_count(args);
                          if (nr_samples < 2) return 0;
                          return sum_of_squares_/(nr_samples-1);
        }

    private:

        Sample previous_mean_;
        Sample sum_of_squares_;
    };
} // namespace impl

///////////////////////////////////////////////////////////////////////////////
// tag:: rolling_variance
//
namespace tag
{
    struct rolling_variance
      : depends_on< rolling_window_plus1, rolling_count, rolling_mean>
    {
        /// INTERNAL ONLY
        ///
        typedef accumulators::impl::rolling_variance_impl< mpl::_1 > impl;

        #ifdef BOOST_ACCUMULATORS_DOXYGEN_INVOKED
        /// tag::rolling_window::window_size named parameter
        static boost::parameter::keyword<tag::rolling_window_size> const window_size;
        #endif
    };
} // namespace tag

///////////////////////////////////////////////////////////////////////////////
// extract::rolling_variance
//
namespace extract
{
    extractor<tag::rolling_variance> const rolling_variance = {};

    BOOST_ACCUMULATORS_IGNORE_GLOBAL(rolling_variance)
}

using extract::rolling_variance;

}} // namespace boost::accumulators

#endif