yat/statistics/Chi2.h

Code
Comments
Other
Rev Date Author Line
3723 05 Feb 18 peter 1 #ifndef _theplu_yat_statistics_chi2_
3723 05 Feb 18 peter 2 #define _theplu_yat_statistics_chi2_
3723 05 Feb 18 peter 3
3723 05 Feb 18 peter 4 // $Id$
3723 05 Feb 18 peter 5
3723 05 Feb 18 peter 6 /*
4207 26 Aug 22 peter 7   Copyright (C) 2018, 2022 Peter Johansson
3723 05 Feb 18 peter 8
3723 05 Feb 18 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
3723 05 Feb 18 peter 10
3723 05 Feb 18 peter 11   The yat library is free software; you can redistribute it and/or
3723 05 Feb 18 peter 12   modify it under the terms of the GNU General Public License as
3723 05 Feb 18 peter 13   published by the Free Software Foundation; either version 3 of the
3723 05 Feb 18 peter 14   License, or (at your option) any later version.
3723 05 Feb 18 peter 15
3723 05 Feb 18 peter 16   The yat library is distributed in the hope that it will be useful,
3723 05 Feb 18 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
3723 05 Feb 18 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3723 05 Feb 18 peter 19   General Public License for more details.
3723 05 Feb 18 peter 20
3723 05 Feb 18 peter 21   You should have received a copy of the GNU General Public License
3723 05 Feb 18 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
3723 05 Feb 18 peter 23 */
3723 05 Feb 18 peter 24
3723 05 Feb 18 peter 25 #include <yat/utility/Matrix.h>
3723 05 Feb 18 peter 26
3723 05 Feb 18 peter 27 namespace theplu {
3723 05 Feb 18 peter 28 namespace yat {
3723 05 Feb 18 peter 29 namespace statistics {
3723 05 Feb 18 peter 30
3723 05 Feb 18 peter 31   /**
3723 05 Feb 18 peter 32      \brief Chi squared test
3723 05 Feb 18 peter 33
3723 05 Feb 18 peter 34      \since New in yat 0.16
3723 05 Feb 18 peter 35   */
3723 05 Feb 18 peter 36   class Chi2
3723 05 Feb 18 peter 37   {
3723 05 Feb 18 peter 38   public:
3723 05 Feb 18 peter 39     /**
3723 05 Feb 18 peter 40        \return degrees of freedom
3723 05 Feb 18 peter 41      */
3723 05 Feb 18 peter 42     unsigned int dof(void) const;
3723 05 Feb 18 peter 43
3723 05 Feb 18 peter 44     /**
3723 05 Feb 18 peter 45        The expected value in element ij is calculated as sum of values
3723 05 Feb 18 peter 46        in row i times sum of values in column j divided by total sum
3723 05 Feb 18 peter 47        of values.
3723 05 Feb 18 peter 48
3723 05 Feb 18 peter 49        \return expected values if there is no association between row
3723 05 Feb 18 peter 50        and column
3723 05 Feb 18 peter 51      */
3723 05 Feb 18 peter 52     const utility::Matrix& expected(void) const;
3723 05 Feb 18 peter 53
3723 05 Feb 18 peter 54     /**
3723 05 Feb 18 peter 55        \return two-sided p-value
3723 05 Feb 18 peter 56     */
3723 05 Feb 18 peter 57     double p_value(void) const;
3723 05 Feb 18 peter 58
3723 05 Feb 18 peter 59     /**
3725 08 Feb 18 peter 60        Randomise the data such that row sums and column sums are
3725 08 Feb 18 peter 61        conserved and calculate chi-squared. Return fraction of time
3725 08 Feb 18 peter 62        the randomised chi-squared is greater (or equal) than the
3725 08 Feb 18 peter 63        chi-squared calculated on original data.
3725 08 Feb 18 peter 64
3725 08 Feb 18 peter 65        \param N number of randomisations
3725 08 Feb 18 peter 66
3725 08 Feb 18 peter 67        Complexity: N x M, where M is total number of samples
3725 08 Feb 18 peter 68      */
3725 08 Feb 18 peter 69     double p_value(size_t N) const;
3725 08 Feb 18 peter 70
3725 08 Feb 18 peter 71     /**
3723 05 Feb 18 peter 72        \return Chi2 sum calculated previously calling operator(void)
3723 05 Feb 18 peter 73      */
3723 05 Feb 18 peter 74     double operator()(void);
3723 05 Feb 18 peter 75
3723 05 Feb 18 peter 76     /**
3723 05 Feb 18 peter 77        Sum is calculated as
3723 05 Feb 18 peter 78
3723 05 Feb 18 peter 79        \f$ sum_{ij} \frac{ (X_{ij}-E_{ij})^2}{E_{ij}} \f$
3723 05 Feb 18 peter 80
3723 05 Feb 18 peter 81        where is X is data and E is matrix with expected values (see
3723 05 Feb 18 peter 82        expected(void));
3723 05 Feb 18 peter 83
3723 05 Feb 18 peter 84        \return Chi2 sum
3723 05 Feb 18 peter 85     */
4125 14 Jan 22 peter 86     double operator()(const utility::MatrixBase& data);
3723 05 Feb 18 peter 87
3723 05 Feb 18 peter 88   private:
3723 05 Feb 18 peter 89     double chi2_;
3723 05 Feb 18 peter 90     utility::Matrix expected_;
3725 08 Feb 18 peter 91     utility::Vector rowsum_;
3725 08 Feb 18 peter 92     utility::Vector colsum_;
3725 08 Feb 18 peter 93     unsigned long int total_;
3725 08 Feb 18 peter 94
4125 14 Jan 22 peter 95     void calculate_expected(const utility::MatrixBase& data);
4125 14 Jan 22 peter 96     void calculate_rowsum(const utility::MatrixBase& data);
4125 14 Jan 22 peter 97     void calculate_colsum(const utility::MatrixBase& data);
3725 08 Feb 18 peter 98
3725 08 Feb 18 peter 99     // Create a matrix with non-negative integers such that the sum of
3725 08 Feb 18 peter 100     // row r is rowsum_(r) and sum of column c is colsum_(c)
3725 08 Feb 18 peter 101     utility::Matrix generate_random_data(void) const;
3725 08 Feb 18 peter 102
4125 14 Jan 22 peter 103     double chi2sum(const utility::MatrixBase& data,
3725 08 Feb 18 peter 104                    const utility::Matrix& expected) const;
3725 08 Feb 18 peter 105
3725 08 Feb 18 peter 106     size_t draw(const utility::Vector& count, unsigned long int N) const;
3725 08 Feb 18 peter 107
3725 08 Feb 18 peter 108     struct Calculator
3725 08 Feb 18 peter 109     {
3725 08 Feb 18 peter 110       double operator()(double x, double m) const;
3725 08 Feb 18 peter 111     };
3723 05 Feb 18 peter 112    };
3723 05 Feb 18 peter 113
3723 05 Feb 18 peter 114 }}} // of namespace statistics, yat, and theplu
3723 05 Feb 18 peter 115
3723 05 Feb 18 peter 116 #endif