averager_base.h

#ifndef _theplu_yat_statistics_average_base
#define _theplu_yat_statistics_average_base

// $Id: averager_base.h 2809 2012-08-06 13:44:47Z peter $

/*
  Copyright (C) 2012 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <boost/concept_check.hpp>

#include <cmath>
#include <limits>

namespace theplu {
namespace yat {
namespace statistics {

  template<class Derived>
  class averager_base
  {
  public:
    averager_base(void) : n_(0), mean_(0)
    {
      BOOST_CONCEPT_ASSERT((boost::DefaultConstructible<Derived>));
    }

    virtual ~averager_base(void)=0;

    void add(double x, long n=1) {static_cast<Derived*>(this)->add_impl(x, n);}

    double mean(void) const
    { return n_ ? mean_ : std::numeric_limits<double>::quiet_NaN(); }

    long n(void) const { return n_; }

    void rescale(double a) { static_cast<Derived*>(this)->rescale_impl(a); }

    void reset(void)
    {
      Derived tmp;
      *static_cast<Derived*>(this) = tmp;
    }

    double sum_x(void)  const { return n_*mean_; }

  protected:
    long int n_;

    double mean_;

    void add1(double x, long int n);
    void add1(double delta);

    double rescale1(double x);
  };


  template<class Derived>
  class averager_base2 : public averager_base<Derived>
  {
  public:
    averager_base2(void) : cm2_(0) {}

    virtual ~averager_base2(void)=0;

    double cv(void) const { return std() / this->mean(); }

    double standard_error(void) const
    { return std::sqrt(variance()/this->n()); }

    double std(void) const { return std::sqrt(variance()); }

    double std(double m) const { return std::sqrt(variance(m)); }

    double sum_xx(void) const
    { return cm2_ + this->n_*this->mean_*this->mean_; }

    double sum_xx_centered(void)  const { return cm2_; }

    double variance(double m) const
    { return sum_xx()/this->n() + m*m - 2*m*this->mean();}

    double variance(void) const { return sum_xx_centered()/this->n_; }

    double variance_unbiased(void) const
    {
      return (this->n_>1) ? sum_xx_centered()/(this->n_-1) :
        std::numeric_limits<double>::quiet_NaN();
    }

  protected:
    double cm2_;

    void add2(double mean, double cm2, long int n);

    void add2(double delta);

    double rescale2(double x);
  private:
  };


  template<class Derived>
  class averager_base3 : public averager_base2<Derived>
  {
  public:
    averager_base3(void) : averager_base2<Derived>(), cm3_(0) {}

    virtual ~averager_base3(void)=0;

    double central_moment3(void) const { return cm3_/this->n_; }

    double skewness(void) const
    { return central_moment3()/std::pow(this->std(),3); }

    double sum_x3_centered(void) const { return cm3_; }

  protected:
    double cm3_;

    void add3(double mean, double cm2, double cm3, long int n);

    void add3(double delta);

    double rescale3(double x);
  };

  template<class Derived>
  class averager_base4 : public averager_base3<Derived>
  {
  public:
    averager_base4(void) : averager_base3<Derived>(), cm4_(0) {}

    virtual ~averager_base4(void)=0;

    double central_moment4(void) const { return cm4_/this->n_; }

    double kurtosis(void) const
    { return central_moment4()/std::pow(this->variance(),2) - 3.0; }

    double sum_x4_centered(void) const { return cm4_; }
  protected:
    double cm4_;

    void add4(double mean, double cm2, double cm3, double cm4, long int n);

    void add4(double delta);

    double rescale4(double x);
  };


  template<class Derived>
  averager_base<Derived>::~averager_base(void) {}
  template<class Derived>
  averager_base2<Derived>::~averager_base2(void) {}
  template<class Derived>
  averager_base3<Derived>::~averager_base3(void) {}
  template<class Derived>
  averager_base4<Derived>::~averager_base4(void) {}

  template<class Derived>
  void averager_base<Derived>::add1(double x, long int n)
  {
    if (n == -this->n_) {
      this->reset();
      return;
    }
    this->mean_ = (this->n_*this->mean_ + n*x) / (this->n_ + n);
    this->n_+=n;
  }


  template<class Derived>
  void averager_base<Derived>::add1(double delta)
  {
    ++this->n_;
    this->mean_ += delta * 1/(this->n_);
  }


  template<class Derived>
  void averager_base2<Derived>::add2(double m, double cm2, long int n)
  {
    if (n == -this->n_) {
      this->reset();
      return;
    }
    const double delta = this->mean_ - m;
    this->cm2_ += cm2 + this->n_*n*delta*delta/(this->n_+n);
    this->add1(m, n);
  }


  template<class Derived>
  void averager_base2<Derived>::add2(double delta)
  {
    this->cm2_ += this->n_*delta*delta/(this->n_+1);
    this->add1(delta);
  }


  template<class Derived>
  void averager_base3<Derived>::add3(double mean, double cm2, double cm3,
                                     long int n)
  {
    if (n == -this->n_) {
      this->reset();
      return;
    }
    const double n_inverse = 1.0/(this->n_+n);
    const double delta = mean - this->mean_;
    const double r = n_inverse*delta;

    this->cm3_ += cm3 + r*r*delta*this->n_*n*(this->n_-n) +
      3*r*(this->n_*cm2 - n*this->cm2_);
    this->add2(mean, cm2, n);
  }


  template<class Derived>
  void averager_base3<Derived>::add3(double delta)
  {
    const double r = delta/(this->n_+1);

    this->cm3_ += r*r*delta*this->n_*(this->n_-1) - 3*r*this->cm2_;
    this->add2(delta);
  }


  template<class Derived>
  void averager_base4<Derived>::add4(double mean, double cm2, double cm3,
                                     double cm4, long int n)
  {
    if (n == -this->n_) {
      this->reset();
      return;
    }
    const double n_inverse = 1.0/(this->n_+n);
    const double delta = mean - this->mean_;
    const double r = n_inverse*delta;

    cm4_ += cm4 +
      std::pow(r,3)*delta*this->n_*n*(this->n_*this->n_-this->n_*n+n*n) +
      6*r*r*(this->n_*this->n_*cm2+n*n*this->cm2_) +
      4*r*(this->n_*cm3-n*this->cm3_);
    this->add3(mean, cm2, cm3, n);
  }


  template<class Derived>
  void averager_base4<Derived>::add4(double delta)
  {
    const double r = delta/(this->n_+1);

    cm4_ +=
      std::pow(r,3)*delta*this->n_*(this->n_*this->n_-this->n_+1) +
      6*r*r*(this->cm2_) +
      4*r*(0-this->cm3_);

    this->add3(delta);
  }


  template<class Derived>
  double averager_base<Derived>::rescale1(double x)
  {
    this->mean_ *= x;
    return x;
  }


  template<class Derived>
  double averager_base2<Derived>::rescale2(double x)
  {
    this->rescale1(x);
    double factor = x*x;
    this->cm2_ *= factor;
    return factor;
  }


  template<class Derived>
  double averager_base3<Derived>::rescale3(double x)
  {
    double factor = x*this->rescale2(x);
    this->cm3_ *= factor;
    return factor;
  }


  template<class Derived>
  double averager_base4<Derived>::rescale4(double x)
  {
    double factor = x*this->rescale3(x);
    this->cm3_ *= factor;
    return factor;
  }

}}}
#endif