test/averager4.cc

Code
Comments
Other
Rev Date Author Line
2799 28 Jul 12 peter 1 // $Id$
2799 28 Jul 12 peter 2
2799 28 Jul 12 peter 3 /*
2799 28 Jul 12 peter 4   Copyright (C) 2012 Peter Johansson
2799 28 Jul 12 peter 5
2799 28 Jul 12 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
2799 28 Jul 12 peter 7
2799 28 Jul 12 peter 8   The yat library is free software; you can redistribute it and/or
2799 28 Jul 12 peter 9   modify it under the terms of the GNU General Public License as
2799 28 Jul 12 peter 10   published by the Free Software Foundation; either version 3 of the
2799 28 Jul 12 peter 11   License, or (at your option) any later version.
2799 28 Jul 12 peter 12
2799 28 Jul 12 peter 13   The yat library is distributed in the hope that it will be useful,
2799 28 Jul 12 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
2799 28 Jul 12 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2799 28 Jul 12 peter 16   General Public License for more details.
2799 28 Jul 12 peter 17
2799 28 Jul 12 peter 18   You should have received a copy of the GNU General Public License
2799 28 Jul 12 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
2799 28 Jul 12 peter 20 */
2799 28 Jul 12 peter 21
2881 18 Nov 12 peter 22 #include <config.h>
2881 18 Nov 12 peter 23
2799 28 Jul 12 peter 24 #include "Suite.h"
2799 28 Jul 12 peter 25
2799 28 Jul 12 peter 26 #include "yat/statistics/Averager4.h"
2799 28 Jul 12 peter 27 #include "yat/statistics/utility.h"
2799 28 Jul 12 peter 28 #include "yat/utility/Vector.h"
2799 28 Jul 12 peter 29
2799 28 Jul 12 peter 30 using namespace theplu::yat;
2799 28 Jul 12 peter 31 using namespace theplu::yat::statistics;
2799 28 Jul 12 peter 32
2799 28 Jul 12 peter 33 using theplu::yat::test::Suite;
2799 28 Jul 12 peter 34
2799 28 Jul 12 peter 35 int main(int argc, char* argv[])
2799 28 Jul 12 peter 36 {
2799 28 Jul 12 peter 37   Suite suite(argc, argv);
2799 28 Jul 12 peter 38
2799 28 Jul 12 peter 39   Averager4 a;
2799 28 Jul 12 peter 40   Averager4 copy(a);
2799 28 Jul 12 peter 41   a = copy;
2799 28 Jul 12 peter 42   a += copy;
2799 28 Jul 12 peter 43
2809 06 Aug 12 peter 44   suite.out() << "testing with n=9\n";
2799 28 Jul 12 peter 45   a.add(0, 9);
2799 28 Jul 12 peter 46   a.add(10);
2799 28 Jul 12 peter 47
2799 28 Jul 12 peter 48   if (!suite.equal(a.central_moment3(), (-9.0 + std::pow(9.0,3))/10.0)) {
2799 28 Jul 12 peter 49     suite.add(false);
2799 28 Jul 12 peter 50     suite.err() << "central moment3 failed\n";
2809 06 Aug 12 peter 51     suite.err() << "sum x3 centered: " << a.sum_x3_centered() << "\n";
2809 06 Aug 12 peter 52     suite.err() << "mean: " << a.mean() << "\n";
2809 06 Aug 12 peter 53     suite.err() << "variance: " << a.variance() << "\n";
2799 28 Jul 12 peter 54   }
2799 28 Jul 12 peter 55   if (!suite.equal(a.central_moment4(), (9 + 9*9*9*9)/10)) {
2799 28 Jul 12 peter 56     suite.add(false);
2799 28 Jul 12 peter 57     suite.err() << "central moment4 failed\n";
2799 28 Jul 12 peter 58   }
2799 28 Jul 12 peter 59   a.reset();
2809 06 Aug 12 peter 60   if (a.n()) {
2809 06 Aug 12 peter 61     suite.add(false);
2809 06 Aug 12 peter 62     suite.err() << "error n: " << a.n() << "\n";
2809 06 Aug 12 peter 63   }
2799 28 Jul 12 peter 64
2809 06 Aug 12 peter 65   suite.out() << "testing with 10 added values\n";
2799 28 Jul 12 peter 66   utility::Vector vec(10);
2799 28 Jul 12 peter 67   for (size_t i=0; i<vec.size(); ++i) {
2803 30 Jul 12 peter 68     vec(i) = i*i;
2799 28 Jul 12 peter 69     a.add(vec(i), 1);
2799 28 Jul 12 peter 70   }
2809 06 Aug 12 peter 71   double m = a.mean();
2809 06 Aug 12 peter 72   suite.out() << "m: " << m << "\n";
2799 28 Jul 12 peter 73
2799 28 Jul 12 peter 74   double sum = 0;
2799 28 Jul 12 peter 75   for (size_t i=0; i<vec.size(); ++i)
2799 28 Jul 12 peter 76     sum += std::pow(vec(i)-m, 3);
2799 28 Jul 12 peter 77   sum = sum/vec.size();
2803 30 Jul 12 peter 78   if (!suite.equal(a.central_moment3(), sum, 2)) {
2799 28 Jul 12 peter 79     suite.add(false);
2799 28 Jul 12 peter 80     suite.err() << "central moment3 failed: expected " << sum << "\n";
2799 28 Jul 12 peter 81   }
2799 28 Jul 12 peter 82
2799 28 Jul 12 peter 83   sum = 0;
2799 28 Jul 12 peter 84   for (size_t i=0; i<vec.size(); ++i)
2799 28 Jul 12 peter 85     sum += std::pow(vec(i)-m, 4)/vec.size();
2799 28 Jul 12 peter 86   if (!suite.equal(a.central_moment4(), sum)) {
2799 28 Jul 12 peter 87     suite.add(false);
2799 28 Jul 12 peter 88     suite.err() << "central moment4 failed\n";
2799 28 Jul 12 peter 89   }
2799 28 Jul 12 peter 90
2799 28 Jul 12 peter 91   suite.out() << "comparing against GSL implementations\n";
2803 30 Jul 12 peter 92   double correct_factor = std::sqrt((vec.size()-1.0)/vec.size());
2799 28 Jul 12 peter 93   if (!suite.add(suite.equal(a.skewness(),
2803 30 Jul 12 peter 94                              skewness(vec)/std::pow(correct_factor,3),
2799 28 Jul 12 peter 95                              10))) {
2799 28 Jul 12 peter 96     suite.err() << "skewness failed\n";
2803 30 Jul 12 peter 97     suite.err() << "Averager4: " << a.skewness() << "\n";
2803 30 Jul 12 peter 98     suite.err() << "GSL:       " << skewness(vec) << "\n";
2803 30 Jul 12 peter 99     suite.err() << "factor: " << std::pow(correct_factor,3)  << "\n";
2799 28 Jul 12 peter 100   }
2799 28 Jul 12 peter 101
2799 28 Jul 12 peter 102   if (!suite.add(suite.equal(a.kurtosis(),
2803 30 Jul 12 peter 103                              (kurtosis(vec)+3)/std::pow(correct_factor,4)-3,
2799 28 Jul 12 peter 104                              10))) {
2799 28 Jul 12 peter 105     suite.err() << "kurtosis failed\n";
2803 30 Jul 12 peter 106     suite.err() << "Averager4: " << a.kurtosis() << "\n";
2803 30 Jul 12 peter 107     suite.err() << "GSL:       " << kurtosis(vec) << "\n";
2803 30 Jul 12 peter 108     suite.err() << "factor: " << std::pow(correct_factor,4)  << "\n";
2799 28 Jul 12 peter 109   }
2799 28 Jul 12 peter 110
2799 28 Jul 12 peter 111   suite.out() << "test operator+=\n";
2799 28 Jul 12 peter 112   Averager4 a1;
2799 28 Jul 12 peter 113   a1.add(0);
2799 28 Jul 12 peter 114   a1.add(-2.5,4);
2799 28 Jul 12 peter 115
2799 28 Jul 12 peter 116   Averager4 a2;
2799 28 Jul 12 peter 117   a2.add(1,4);
2799 28 Jul 12 peter 118   a2.add(6,1);
2799 28 Jul 12 peter 119
2809 06 Aug 12 peter 120   suite.out() << a2.mean() << "\n";
2803 30 Jul 12 peter 121   // cm3 = (-4*1 + 4^3)/5 = 60/5 = 12
2803 30 Jul 12 peter 122   if (!suite.add(suite.equal(a2.central_moment3(), 12)))
2799 28 Jul 12 peter 123     suite.err() << "central moment3 failed\n";
2799 28 Jul 12 peter 124
2799 28 Jul 12 peter 125   a1 += a2;
2809 06 Aug 12 peter 126   if (!suite.add(suite.equal(a1.mean(), 0)))
2799 28 Jul 12 peter 127     suite.err() << "mean failed\n";
2809 06 Aug 12 peter 128   if (!suite.add(suite.equal(a1.variance(),
2799 28 Jul 12 peter 129                              (4*2.5*2.5+4*1*1+1*6*6)/10.0)))
2799 28 Jul 12 peter 130     suite.err() << "variance failed\n";
2799 28 Jul 12 peter 131   if (!suite.add(suite.equal(a1.central_moment3(),
2803 30 Jul 12 peter 132                              (4*std::pow(-2.5,3)+4+std::pow(6.0,3))/10)))
2799 28 Jul 12 peter 133     suite.err() << "central_moment3 failed\n";
2799 28 Jul 12 peter 134   if (!suite.add(suite.equal(a1.central_moment4(),
2799 28 Jul 12 peter 135                              (4*std::pow(2.5,4)+4+std::pow(6.0,4))/10)))
2799 28 Jul 12 peter 136     suite.err() << "central_moment4 failed\n";
2799 28 Jul 12 peter 137
2809 06 Aug 12 peter 138   a2.rescale(2);
2799 28 Jul 12 peter 139
2799 28 Jul 12 peter 140   return suite.return_value();
2799 28 Jul 12 peter 141 }

Generated on Sat Oct 14 02:20:38 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)