source: trunk/test/averager4.cc @ 2803

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

fixed spurious test failures as well as a bug in operator+=. refs #705

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.8 KB
Line 
1// $Id: averager4.cc 2803 2012-07-30 05:18:25Z 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 "Suite.h"
23
24#include "yat/statistics/Averager4.h"
25#include "yat/statistics/utility.h"
26#include "yat/utility/Vector.h"
27
28using namespace theplu::yat;
29using namespace theplu::yat::statistics;
30
31using theplu::yat::test::Suite;
32
33int main(int argc, char* argv[])
34{
35  Suite suite(argc, argv);
36
37  Averager4 a;
38  Averager4 copy(a);
39  a = copy;
40  a += copy;
41
42  a.add(0, 9);
43  a.add(10);
44
45  if (!suite.equal(a.central_moment3(), (-9.0 + std::pow(9.0,3))/10.0)) {
46    suite.add(false);
47    suite.err() << "central moment3 failed\n";
48  }
49  if (!suite.equal(a.central_moment4(), (9 + 9*9*9*9)/10)) {
50    suite.add(false);
51    suite.err() << "central moment4 failed\n";
52  }
53  a.reset();
54
55  suite.out() << "=========================\n";
56  utility::Vector vec(10);
57  for (size_t i=0; i<vec.size(); ++i) {
58    vec(i) = i*i;
59    a.add(vec(i), 1);
60  }
61  double m = a.averager().mean();
62
63  double sum = 0;
64  for (size_t i=0; i<vec.size(); ++i)
65    sum += std::pow(vec(i)-m, 3);
66  sum = sum/vec.size();
67  if (!suite.equal(a.central_moment3(), sum, 2)) {
68    suite.add(false);
69    suite.err() << "central moment3 failed: expected " << sum << "\n";
70  }
71  //  return suite.return_value();
72
73  sum = 0;
74  for (size_t i=0; i<vec.size(); ++i)
75    sum += std::pow(vec(i)-m, 4)/vec.size();
76  if (!suite.equal(a.central_moment4(), sum)) {
77    suite.add(false);
78    suite.err() << "central moment4 failed\n";
79  }
80
81  suite.out() << "comparing against GSL implementations\n";
82  double correct_factor = std::sqrt((vec.size()-1.0)/vec.size());
83  if (!suite.add(suite.equal(a.skewness(),
84                             skewness(vec)/std::pow(correct_factor,3),
85                             10))) {
86    suite.err() << "skewness failed\n";
87    suite.err() << "Averager4: " << a.skewness() << "\n";
88    suite.err() << "GSL:       " << skewness(vec) << "\n";
89    suite.err() << "factor: " << std::pow(correct_factor,3) << "\n";
90  }
91
92  if (!suite.add(suite.equal(a.kurtosis(),
93                             (kurtosis(vec)+3)/std::pow(correct_factor,4)-3,
94                             10))) {
95    suite.err() << "kurtosis failed\n";
96    suite.err() << "Averager4: " << a.kurtosis() << "\n";
97    suite.err() << "GSL:       " << kurtosis(vec) << "\n";
98    suite.err() << "factor: " << std::pow(correct_factor,4) << "\n";
99
100    double var = a.averager().variance_unbiased();
101    suite.err() << a.central_moment4()/(var*var) - 3 << "\n";
102  }
103
104  suite.out() << "test operator+=\n";
105  Averager4 a1;
106  a1.add(0);
107  a1.add(-2.5,4);
108
109  Averager4 a2;
110  a2.add(1,4);
111  a2.add(6,1);
112
113  // cm3 = (-4*1 + 4^3)/5 = 60/5 = 12
114  if (!suite.add(suite.equal(a2.central_moment3(), 12)))
115    suite.err() << "central moment3 failed\n";
116
117  a1 += a2;
118  if (!suite.add(suite.equal(a1.averager().mean(), 0)))
119    suite.err() << "mean failed\n";
120  if (!suite.add(suite.equal(a1.averager().variance(),
121                             (4*2.5*2.5+4*1*1+1*6*6)/10.0)))
122    suite.err() << "variance failed\n";
123  if (!suite.add(suite.equal(a1.central_moment3(),
124                             (4*std::pow(-2.5,3)+4+std::pow(6.0,3))/10)))
125    suite.err() << "central_moment3 failed\n";
126  if (!suite.add(suite.equal(a1.central_moment4(),
127                             (4*std::pow(2.5,4)+4+std::pow(6.0,4))/10)))
128    suite.err() << "central_moment4 failed\n";
129
130
131  Averager aver = a.averager();
132  aver.mean();
133
134  return suite.return_value();
135}
Note: See TracBrowser for help on using the repository browser.