source: trunk/test/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 "Suite.h"

#include "yat/statistics/Averager4.h"
#include "yat/statistics/utility.h"
#include "yat/utility/Vector.h"

using namespace theplu::yat;
using namespace theplu::yat::statistics;

using theplu::yat::test::Suite;

int main(int argc, char* argv[])
{
  Suite suite(argc, argv);

  Averager4 a;
  Averager4 copy(a);
  a = copy;
  a += copy;

  a.add(0, 9);
  a.add(10);

  if (!suite.equal(a.central_moment3(), (-9.0 + std::pow(9.0,3))/10.0)) {
    suite.add(false);
    suite.err() << "central moment3 failed\n";
  }
  if (!suite.equal(a.central_moment4(), (9 + 9*9*9*9)/10)) {
    suite.add(false);
    suite.err() << "central moment4 failed\n";
  }
  a.reset();

  suite.out() << "=========================\n";
  utility::Vector vec(10);
  for (size_t i=0; i<vec.size(); ++i) {
    vec(i) = i*i;
    a.add(vec(i), 1);
  }
  double m = a.averager().mean();

  double sum = 0;
  for (size_t i=0; i<vec.size(); ++i)
    sum += std::pow(vec(i)-m, 3);
  sum = sum/vec.size();
  if (!suite.equal(a.central_moment3(), sum, 2)) {
    suite.add(false);
    suite.err() << "central moment3 failed: expected " << sum << "\n";
  }
  //  return suite.return_value();

  sum = 0;
  for (size_t i=0; i<vec.size(); ++i)
    sum += std::pow(vec(i)-m, 4)/vec.size();
  if (!suite.equal(a.central_moment4(), sum)) {
    suite.add(false);
    suite.err() << "central moment4 failed\n";
  }

  suite.out() << "comparing against GSL implementations\n";
  double correct_factor = std::sqrt((vec.size()-1.0)/vec.size());
  if (!suite.add(suite.equal(a.skewness(),
                             skewness(vec)/std::pow(correct_factor,3),
                             10))) {
    suite.err() << "skewness failed\n";
    suite.err() << "Averager4: " << a.skewness() << "\n";
    suite.err() << "GSL:       " << skewness(vec) << "\n";
    suite.err() << "factor: " << std::pow(correct_factor,3) << "\n";
  }

  if (!suite.add(suite.equal(a.kurtosis(),
                             (kurtosis(vec)+3)/std::pow(correct_factor,4)-3,
                             10))) {
    suite.err() << "kurtosis failed\n";
    suite.err() << "Averager4: " << a.kurtosis() << "\n";
    suite.err() << "GSL:       " << kurtosis(vec) << "\n";
    suite.err() << "factor: " << std::pow(correct_factor,4) << "\n";

    double var = a.averager().variance_unbiased();
    suite.err() << a.central_moment4()/(var*var) - 3 << "\n";
  }

  suite.out() << "test operator+=\n";
  Averager4 a1;
  a1.add(0);
  a1.add(-2.5,4);

  Averager4 a2;
  a2.add(1,4);
  a2.add(6,1);

  // cm3 = (-4*1 + 4^3)/5 = 60/5 = 12
  if (!suite.add(suite.equal(a2.central_moment3(), 12)))
    suite.err() << "central moment3 failed\n";

  a1 += a2;
  if (!suite.add(suite.equal(a1.averager().mean(), 0)))
    suite.err() << "mean failed\n";
  if (!suite.add(suite.equal(a1.averager().variance(),
                             (4*2.5*2.5+4*1*1+1*6*6)/10.0)))
    suite.err() << "variance failed\n";
  if (!suite.add(suite.equal(a1.central_moment3(),
                             (4*std::pow(-2.5,3)+4+std::pow(6.0,3))/10)))
    suite.err() << "central_moment3 failed\n";
  if (!suite.add(suite.equal(a1.central_moment4(),
                             (4*std::pow(2.5,4)+4+std::pow(6.0,4))/10)))
    suite.err() << "central_moment4 failed\n";


  Averager aver = a.averager();
  aver.mean();

  return suite.return_value();
}