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

// $Id: Anova.cc 4114 2021-10-13 04:36:16Z peter $

/*
  Copyright (C) 2021 Peter Johansson

  This file is part of the yat library, https://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <https://www.gnu.org/licenses/>.
*/

#include "Anova.h"

#include <gsl/gsl_cdf.h>

#include <cassert>

namespace theplu {
namespace yat {
namespace statistics {

  Anova::Anova(size_t n)
    : aver_(n)
  {
  }


  void Anova::add(double x, size_t idx, long int n)
  {
    assert(idx < aver_.size());
    aver_[idx].add(x, n);
    total_.add(x, n);
  }


  double Anova::F(void) const
  {
    double intra_ss = 0;
    for (const auto& a : aver_)
      intra_ss += a.sum_xx_centered();
    double total_ss = total_.sum_xx_centered();

    /*
      Remember that
      ss = \sum (x-m)^2 = \sum (x-m_i+m_i-m)^2 =
      \sum [(x-m_i)^2 + (m_i-m)^2] =
      \sum [ss_i + n_i (m_i-m)^2]
      where the last sum runs over groups and the other over samples, so
      total ss = intra ss + inter ss
     */
    double inter_ss = total_ss - intra_ss;

    // return (inter_ss / inter_df) / (intra_ss / intra_df)
    return (inter_ss * intra_df()) / (intra_ss * inter_df());
  }


  size_t Anova::inter_df(void) const
  {
    return aver_.size() - 1;
  }


  size_t Anova::intra_df(void) const
  {
    return total_.n() - aver_.size();
  }


  double Anova::p_value(void) const
  {
    return gsl_cdf_fdist_Q(F(), inter_df(), intra_df());
  }


  void Anova::reset(void)
  {
    for (auto& a : aver_)
      a.reset();
    total_.reset();
  }

}}}