source: trunk/yat/statistics/Anova.cc @ 4114

Last change on this file since 4114 was 4114, checked in by Peter, 9 months ago

new class for one-way ANOVA test

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 1.9 KB
Line 
1// $Id: Anova.cc 4114 2021-10-13 04:36:16Z peter $
2
3/*
4  Copyright (C) 2021 Peter Johansson
5
6  This file is part of the yat library, https://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 <https://www.gnu.org/licenses/>.
20*/
21
22#include "Anova.h"
23
24#include <gsl/gsl_cdf.h>
25
26#include <cassert>
27
28namespace theplu {
29namespace yat {
30namespace statistics {
31
32  Anova::Anova(size_t n)
33    : aver_(n)
34  {
35  }
36
37
38  void Anova::add(double x, size_t idx, long int n)
39  {
40    assert(idx < aver_.size());
41    aver_[idx].add(x, n);
42    total_.add(x, n);
43  }
44
45
46  double Anova::F(void) const
47  {
48    double intra_ss = 0;
49    for (const auto& a : aver_)
50      intra_ss += a.sum_xx_centered();
51    double total_ss = total_.sum_xx_centered();
52
53    /*
54      Remember that
55      ss = \sum (x-m)^2 = \sum (x-m_i+m_i-m)^2 =
56      \sum [(x-m_i)^2 + (m_i-m)^2] =
57      \sum [ss_i + n_i (m_i-m)^2]
58      where the last sum runs over groups and the other over samples, so
59      total ss = intra ss + inter ss
60     */
61    double inter_ss = total_ss - intra_ss;
62
63    // return (inter_ss / inter_df) / (intra_ss / intra_df)
64    return (inter_ss * intra_df()) / (intra_ss * inter_df());
65  }
66
67
68  size_t Anova::inter_df(void) const
69  {
70    return aver_.size() - 1;
71  }
72
73
74  size_t Anova::intra_df(void) const
75  {
76    return total_.n() - aver_.size();
77  }
78
79
80  double Anova::p_value(void) const
81  {
82    return gsl_cdf_fdist_Q(F(), inter_df(), intra_df());
83  }
84
85
86  void Anova::reset(void)
87  {
88    for (auto& a : aver_)
89      a.reset();
90    total_.reset();
91  }
92
93}}}
Note: See TracBrowser for help on using the repository browser.