source: trunk/yat/statistics/Averager4.cc @ 2803

// $Id: Averager4.cc 2803 2012-07-30 05:18:25Z 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 "Averager4.h"

#include <cassert>
#include <cmath>

namespace theplu {
namespace yat {
namespace statistics {

  Averager4::Averager4(void)
    : cm3_(0), cm4_(0) {}


  void Averager4::add(double x, long int n)
  {
    add(x, 0, 0, 0, n);
    averager_.add(x, n);
  }


  void Averager4::add(double mean, double cm2, double cm3, double cm4,
                      long int n2)
  {
    // rhs empty - nothing to do
    if (!n2)
      return;
    const double n1 = averager_.n();
    // lhs empty - simple assignment
    if (!n1) {
      cm3_ = cm3;
      cm4_ = cm4;
      return;
    }
    const double n_inverse = 1.0/(n1+n2);
    const double delta = mean - averager_.mean();
    const double r = n_inverse*delta;

    double new_cm3 = cm3_ + cm3 + r*r*delta*n1*n2*(n1-n2) +
      3*r*(n1*cm2 - n2*averager_.sum_xx_centered());

    cm4_ += cm4 + std::pow(r,3)*delta*n1*n2*(n1*n1-n1*n2+n2*n2) +
      6*r*r*(n1*n1*cm2+n2*n2*averager().sum_xx_centered()) +
      4*r*(n1*cm3-n2*cm3_);

    cm3_ = new_cm3;
  }


  const Averager& Averager4::averager(void) const
  {
    return averager_;
  }


  double Averager4::central_moment3(void) const
  {
    return cm3_/averager_.n();
  }


  double Averager4::central_moment4(void) const
  {
    return cm4_/averager().n();
  }


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


  void Averager4::reset(void)
  {
    averager_.reset();
    cm3_ = 0;
    cm4_ = 0;
  }


  double Averager4::skewness(void) const
  {
    return central_moment3() / std::pow(averager_.std(), 3);
  }


  Averager4& Averager4::operator+=(const Averager4& rhs)
  {
    add(rhs.averager().mean(), rhs.averager().sum_xx_centered(),
        rhs.cm3_, rhs.cm4_, rhs.averager().n());
    averager_ += rhs.averager_;
    return *this;
  }

}}} // of namespace statistics, yat, and theplu