yat/statistics/LikelihoodRatioTestBinomial.cc

Code
Comments
Other
Rev Date Author Line
3134 27 Nov 13 peter 1 // $Id$
3134 27 Nov 13 peter 2
3134 27 Nov 13 peter 3 /*
3134 27 Nov 13 peter 4   Copyright (C) 2013 Peter Johansson
3134 27 Nov 13 peter 5
3134 27 Nov 13 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3134 27 Nov 13 peter 7
3134 27 Nov 13 peter 8   The yat library is free software; you can redistribute it and/or
3134 27 Nov 13 peter 9   modify it under the terms of the GNU General Public License as
3134 27 Nov 13 peter 10   published by the Free Software Foundation; either version 3 of the
3134 27 Nov 13 peter 11   License, or (at your option) any later version.
3134 27 Nov 13 peter 12
3134 27 Nov 13 peter 13   The yat library is distributed in the hope that it will be useful,
3134 27 Nov 13 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3134 27 Nov 13 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3134 27 Nov 13 peter 16   General Public License for more details.
3134 27 Nov 13 peter 17
3134 27 Nov 13 peter 18   You should have received a copy of the GNU General Public License
3134 27 Nov 13 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3134 27 Nov 13 peter 20 */
3134 27 Nov 13 peter 21
3134 27 Nov 13 peter 22 #include <config.h>
3134 27 Nov 13 peter 23
3138 29 Nov 13 peter 24 #include "LikelihoodRatioTestBinomial.h"
3134 27 Nov 13 peter 25
3138 29 Nov 13 peter 26 #include <cassert>
3138 29 Nov 13 peter 27 #include <cmath>
3134 27 Nov 13 peter 28 #include <numeric>
3134 27 Nov 13 peter 29
3134 27 Nov 13 peter 30 namespace theplu {
3134 27 Nov 13 peter 31 namespace yat {
3134 27 Nov 13 peter 32 namespace statistics {
3134 27 Nov 13 peter 33
3138 29 Nov 13 peter 34   LikelihoodRatioTestBinomial::LikelihoodRatioTestBinomial(size_t n)
3138 29 Nov 13 peter 35     : count_(2, n, 0)
3134 27 Nov 13 peter 36   {}
3134 27 Nov 13 peter 37
3138 29 Nov 13 peter 38   void LikelihoodRatioTestBinomial::add(size_t population, bool positive,
3134 27 Nov 13 peter 39                                            long int n)
3134 27 Nov 13 peter 40   {
3138 29 Nov 13 peter 41     assert(population < count_.columns());
3138 29 Nov 13 peter 42     assert(population < count_.columns());
3138 29 Nov 13 peter 43     count_(positive ? 1 : 0, population) += n;
3134 27 Nov 13 peter 44   }
3134 27 Nov 13 peter 45
3134 27 Nov 13 peter 46
3138 29 Nov 13 peter 47   unsigned long int LikelihoodRatioTestBinomial::n(void) const
3134 27 Nov 13 peter 48   {
3138 29 Nov 13 peter 49     return positive() + negative();
3134 27 Nov 13 peter 50   }
3134 27 Nov 13 peter 51
3134 27 Nov 13 peter 52
3138 29 Nov 13 peter 53   unsigned long int LikelihoodRatioTestBinomial::n(size_t i) const
3134 27 Nov 13 peter 54   {
3138 29 Nov 13 peter 55     return positive(i) + negative(i);
3134 27 Nov 13 peter 56   }
3134 27 Nov 13 peter 57
3134 27 Nov 13 peter 58
3138 29 Nov 13 peter 59   unsigned long int LikelihoodRatioTestBinomial::negative(void) const
3134 27 Nov 13 peter 60   {
3138 29 Nov 13 peter 61     return sum(count_.row_const_view(0));
3134 27 Nov 13 peter 62   }
3134 27 Nov 13 peter 63
3134 27 Nov 13 peter 64
3138 29 Nov 13 peter 65   unsigned long int LikelihoodRatioTestBinomial::negative(size_t i) const
3134 27 Nov 13 peter 66   {
3138 29 Nov 13 peter 67     assert(i<count_.columns());
3138 29 Nov 13 peter 68     return count_(0, i);
3134 27 Nov 13 peter 69   }
3134 27 Nov 13 peter 70
3134 27 Nov 13 peter 71
3138 29 Nov 13 peter 72   unsigned long int LikelihoodRatioTestBinomial::positive(void) const
3134 27 Nov 13 peter 73   {
3138 29 Nov 13 peter 74     return sum(count_.row_const_view(1));
3134 27 Nov 13 peter 75   }
3134 27 Nov 13 peter 76
3134 27 Nov 13 peter 77
3138 29 Nov 13 peter 78   unsigned long int LikelihoodRatioTestBinomial::positive(size_t i) const
3134 27 Nov 13 peter 79   {
3138 29 Nov 13 peter 80     assert(i<count_.columns());
3138 29 Nov 13 peter 81     return count_(1, i);
3134 27 Nov 13 peter 82   }
3134 27 Nov 13 peter 83
3134 27 Nov 13 peter 84
3138 29 Nov 13 peter 85   size_t LikelihoodRatioTestBinomial::size(void) const
3134 27 Nov 13 peter 86   {
3138 29 Nov 13 peter 87     return count_.columns();
3134 27 Nov 13 peter 88   }
3134 27 Nov 13 peter 89
3134 27 Nov 13 peter 90
3138 29 Nov 13 peter 91   double LikelihoodRatioTestBinomial::llr(void) const
3134 27 Nov 13 peter 92   {
3138 29 Nov 13 peter 93     double res = 0;
3138 29 Nov 13 peter 94     for (size_t i=0; i<size(); ++i) {
3138 29 Nov 13 peter 95       double ni = n(i);
3138 29 Nov 13 peter 96       if (positive(i))
3138 29 Nov 13 peter 97         res += positive(i) * std::log(positive(i)/ni);
3138 29 Nov 13 peter 98       if (negative(i))
3138 29 Nov 13 peter 99         res += negative(i) * std::log(negative(i)/ni);
3138 29 Nov 13 peter 100     }
3138 29 Nov 13 peter 101     double npos = positive();
3138 29 Nov 13 peter 102     double nneg = negative();
3138 29 Nov 13 peter 103     double n = npos + nneg;
3138 29 Nov 13 peter 104     res -= npos * std::log(npos/n);
3138 29 Nov 13 peter 105     res -= nneg * std::log(nneg/n);
3138 29 Nov 13 peter 106     return res;
3134 27 Nov 13 peter 107   }
3134 27 Nov 13 peter 108
3134 27 Nov 13 peter 109 }}} // of namespace theplu yat statistics