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

Last change on this file since 2799 was 2799, checked in by Peter, 11 years ago

first version of class Averager4. refs #705

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 2.4 KB
Line 
1// $Id: Averager4.cc 2799 2012-07-28 08:17:15Z peter $
2
3/*
4  Copyright (C) 2012 Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with yat. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include "Averager4.h"
23
24#include <cassert>
25#include <cmath>
26
27namespace theplu {
28namespace yat {
29namespace statistics {
30
31  Averager4::Averager4(void)
32    : cm3_(0), cm4_(0) {}
33
34
35  void Averager4::add(double x, long int n)
36  {
37    add(x, 0, 0, 0, n);
38    averager_.add(x, n);
39  }
40
41
42  void Averager4::add(double mean, double cm2, double cm3, double cm4,
43                      long int n2)
44  {
45    // rhs empty - nothing to do
46    if (!n2)
47      return;
48    const double n1 = averager_.n();
49    // lhs empty - simple assignment
50    if (!n1) {
51      cm3_ = cm3;
52      cm4_ = cm4;
53      return;
54    }
55    const double n_inverse = 1.0/(n1+n2);
56    const double delta = mean - averager_.mean();
57    const double r = n_inverse*delta;
58
59    double new_cm3 = cm3_ + cm3 + r*r*delta*n1*n2*(n1-n2) +
60      3*r*(n1*cm2 - n2*averager_.sum_xx_centered());
61
62    cm4_ += cm4 + std::pow(r,3)*delta*n1*n2*(n1*n1-n1*n2+n2*n2) +
63      6*r*r*(n1*n1*cm2+n2*n2*averager().sum_xx_centered()) +
64      4*r*(n1*cm3-n2*cm3_);
65
66    cm3_ = new_cm3;
67  }
68
69
70  const Averager& Averager4::averager(void) const
71  {
72    return averager_;
73  }
74
75
76  double Averager4::central_moment3(void) const
77  {
78    return cm3_/averager_.n();
79  }
80
81
82  double Averager4::central_moment4(void) const
83  {
84    return cm4_/averager().n();
85  }
86
87
88  double Averager4::kurtosis(void) const
89  {
90    return central_moment4() / std::pow(averager_.variance(), 2) - 3.0;
91  }
92
93
94  void Averager4::reset(void)
95  {
96    averager_.reset();
97    cm3_ = 0;
98    cm4_ = 0;
99  }
100
101
102  double Averager4::skewness(void) const
103  {
104    return central_moment3() / std::pow(averager_.std(), 3);
105  }
106
107
108  Averager4& Averager4::operator+=(const Averager4& rhs)
109  {
110    add(rhs.averager().mean(), rhs.averager().variance(),
111        rhs.cm3_, rhs.cm4_, rhs.averager().n());
112    averager_ += rhs.averager_;
113    return *this;
114  }
115
116}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.