yat/statistics/Anova.cc

Code
Comments
Other
Rev Date Author Line
4114 13 Oct 21 peter 1 // $Id$
4114 13 Oct 21 peter 2
4114 13 Oct 21 peter 3 /*
4114 13 Oct 21 peter 4   Copyright (C) 2021 Peter Johansson
4114 13 Oct 21 peter 5
4114 13 Oct 21 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4114 13 Oct 21 peter 7
4114 13 Oct 21 peter 8   The yat library is free software; you can redistribute it and/or
4114 13 Oct 21 peter 9   modify it under the terms of the GNU General Public License as
4114 13 Oct 21 peter 10   published by the Free Software Foundation; either version 3 of the
4114 13 Oct 21 peter 11   License, or (at your option) any later version.
4114 13 Oct 21 peter 12
4114 13 Oct 21 peter 13   The yat library is distributed in the hope that it will be useful,
4114 13 Oct 21 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4114 13 Oct 21 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4114 13 Oct 21 peter 16   General Public License for more details.
4114 13 Oct 21 peter 17
4114 13 Oct 21 peter 18   You should have received a copy of the GNU General Public License
4114 13 Oct 21 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4114 13 Oct 21 peter 20 */
4114 13 Oct 21 peter 21
4115 13 Oct 21 peter 22 #include <config.h>
4115 13 Oct 21 peter 23
4114 13 Oct 21 peter 24 #include "Anova.h"
4114 13 Oct 21 peter 25
4114 13 Oct 21 peter 26 #include <gsl/gsl_cdf.h>
4114 13 Oct 21 peter 27
4114 13 Oct 21 peter 28 #include <cassert>
4114 13 Oct 21 peter 29
4114 13 Oct 21 peter 30 namespace theplu {
4114 13 Oct 21 peter 31 namespace yat {
4114 13 Oct 21 peter 32 namespace statistics {
4114 13 Oct 21 peter 33
4114 13 Oct 21 peter 34   Anova::Anova(size_t n)
4114 13 Oct 21 peter 35     : aver_(n)
4114 13 Oct 21 peter 36   {
4114 13 Oct 21 peter 37   }
4114 13 Oct 21 peter 38
4114 13 Oct 21 peter 39
4114 13 Oct 21 peter 40   void Anova::add(double x, size_t idx, long int n)
4114 13 Oct 21 peter 41   {
4114 13 Oct 21 peter 42     assert(idx < aver_.size());
4114 13 Oct 21 peter 43     aver_[idx].add(x, n);
4114 13 Oct 21 peter 44     total_.add(x, n);
4114 13 Oct 21 peter 45   }
4114 13 Oct 21 peter 46
4114 13 Oct 21 peter 47
4114 13 Oct 21 peter 48   double Anova::F(void) const
4114 13 Oct 21 peter 49   {
4114 13 Oct 21 peter 50     double intra_ss = 0;
4114 13 Oct 21 peter 51     for (const auto& a : aver_)
4114 13 Oct 21 peter 52       intra_ss += a.sum_xx_centered();
4114 13 Oct 21 peter 53     double total_ss = total_.sum_xx_centered();
4114 13 Oct 21 peter 54
4114 13 Oct 21 peter 55     /*
4114 13 Oct 21 peter 56       Remember that
4114 13 Oct 21 peter 57       ss = \sum (x-m)^2 = \sum (x-m_i+m_i-m)^2 =
4114 13 Oct 21 peter 58       \sum [(x-m_i)^2 + (m_i-m)^2] =
4114 13 Oct 21 peter 59       \sum [ss_i + n_i (m_i-m)^2]
4114 13 Oct 21 peter 60       where the last sum runs over groups and the other over samples, so
4114 13 Oct 21 peter 61       total ss = intra ss + inter ss
4114 13 Oct 21 peter 62      */
4114 13 Oct 21 peter 63     double inter_ss = total_ss - intra_ss;
4114 13 Oct 21 peter 64
4114 13 Oct 21 peter 65     // return (inter_ss / inter_df) / (intra_ss / intra_df)
4114 13 Oct 21 peter 66     return (inter_ss * intra_df()) / (intra_ss * inter_df());
4114 13 Oct 21 peter 67   }
4114 13 Oct 21 peter 68
4114 13 Oct 21 peter 69
4114 13 Oct 21 peter 70   size_t Anova::inter_df(void) const
4114 13 Oct 21 peter 71   {
4114 13 Oct 21 peter 72     return aver_.size() - 1;
4114 13 Oct 21 peter 73   }
4114 13 Oct 21 peter 74
4114 13 Oct 21 peter 75
4114 13 Oct 21 peter 76   size_t Anova::intra_df(void) const
4114 13 Oct 21 peter 77   {
4114 13 Oct 21 peter 78     return total_.n() - aver_.size();
4114 13 Oct 21 peter 79   }
4114 13 Oct 21 peter 80
4114 13 Oct 21 peter 81
4114 13 Oct 21 peter 82   double Anova::p_value(void) const
4114 13 Oct 21 peter 83   {
4114 13 Oct 21 peter 84     return gsl_cdf_fdist_Q(F(), inter_df(), intra_df());
4114 13 Oct 21 peter 85   }
4114 13 Oct 21 peter 86
4114 13 Oct 21 peter 87
4114 13 Oct 21 peter 88   void Anova::reset(void)
4114 13 Oct 21 peter 89   {
4114 13 Oct 21 peter 90     for (auto& a : aver_)
4114 13 Oct 21 peter 91       a.reset();
4114 13 Oct 21 peter 92     total_.reset();
4114 13 Oct 21 peter 93   }
4114 13 Oct 21 peter 94
4114 13 Oct 21 peter 95 }}}

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