source: trunk/yat/statistics/Fisher.cc @ 676

Last change on this file since 676 was 675, checked in by Jari Häkkinen, 15 years ago

References #83. Changing project name to yat. Compilation will fail in this revision.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.5 KB
Line 
1// $Id: Fisher.cc 675 2006-10-10 12:08:45Z jari $
2
3/*
4  Copyright (C) The authors contributing to this file.
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/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 2 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 this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "yat/statistics/Fisher.h"
25#include "yat/statistics/Score.h"
26#include "yat/statistics/utility.h"
27#include "yat/classifier/DataLookupWeighted1D.h"
28#include "yat/classifier/Target.h"
29
30#include <gsl/gsl_cdf.h>
31#include <gsl/gsl_randist.h>
32
33namespace theplu {
34namespace statistics {
35
36 
37  Fisher::Fisher(bool b)
38    : Score(b), a_(0), b_(0), c_(0), d_(0), oddsratio_(1.0), value_cutoff_(0)
39  {
40  }
41
42
43  double Fisher::Chi2() const
44  {
45    double a,b,c,d;
46    expected(a,b,c,d);
47    return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b + 
48      (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
49  }
50
51
52  void Fisher::expected(double& a, double& b, double& c, double& d) const
53  {
54    double N = a_+b_+c_+d_;
55    a =((a_+b_)*(a_+c_)) / N;
56    b =((a_+b_)*(b_+d_)) / N;
57    c =((c_+d_)*(a_+c_)) / N;
58    d =((c_+d_)*(b_+d_)) / N;
59  }
60
61  double Fisher::oddsratio(const double a,
62                           const double b,
63                           const double c,
64                           const double d) 
65  {
66    // If a column sum or a row sum is zero, the table is nonsense
67    if ((a==0 || d==0) && (c==0 || b==0)){
68      //Peter, should throw exception
69      std::cerr << "Warning: Fisher: Table is not valid\n";
70      return oddsratio_ = 1.0;
71    }
72
73    oddsratio_=(a*d)/(b*d);
74    if (absolute_ && oddsratio_<1)
75      return 1/oddsratio_;
76    return oddsratio_;
77  }
78
79
80  double Fisher::p_value() const
81  {
82    double p=1;
83    if ( ( a_<minimum_size_ || b_<minimum_size_ || 
84           c_<minimum_size_ || d_<minimum_size_) && !weighted_ ){
85
86      p=p_value_exact();
87    }
88    else
89      p=p_value_approximative();
90
91    if (!absolute_){
92      p=p/2;
93      if (oddsratio_<0.5){
94        // last term because >= not equal to !(<=)
95        u_int a = static_cast<u_int>(a_);
96        u_int b = static_cast<u_int>(b_);
97        u_int c = static_cast<u_int>(c_);
98        u_int d = static_cast<u_int>(d_);
99        return 1-p+gsl_ran_hypergeometric_pdf(a, a+b, c+d, a+c);
100      }
101    }
102    return p;
103  }
104
105
106  double Fisher::p_value_approximative() const
107  {
108    return gsl_cdf_chisq_Q(Chi2(), 1.0);
109  }
110
111  double Fisher::p_value_exact() const
112  { 
113    u_int a = static_cast<u_int>(a_);
114    u_int b = static_cast<u_int>(b_);
115    u_int c = static_cast<u_int>(c_);
116    u_int d = static_cast<u_int>(d_);
117   
118    // Since the calculation is symmetric and cdf_hypergeometric_P
119    // loops up to k we choose the smallest number to be k and mirror
120    // the matrix. This choice makes the p-value two-sided.
121
122    if (a<b && a<c && a<d)
123      return statistics::cdf_hypergeometric_P(a,a+b,c+d,a+c);
124    else if (b<a && b<c && b<d)
125      return  statistics::cdf_hypergeometric_P(b,a+b,c+d,b+d);
126    else if (c<a && c<b && c<d)
127      return  statistics::cdf_hypergeometric_P(c,c+d,a+b,a+c);
128
129    return statistics::cdf_hypergeometric_P(d,c+d,a+b,b+d);
130
131  }
132
133
134    double Fisher::score(const classifier::Target& target, 
135                         const utility::vector& value)
136  {
137    weighted_=false;
138    a_=b_=c_=d_=0;
139    for (size_t i=0; i<target.size(); i++)
140      if (target.binary(i))
141        if (value(i)>value_cutoff_)
142          a_++;
143        else
144          c_++;
145      else
146        if (value(i)>value_cutoff_)
147          b_++;
148        else
149          d_++;
150       
151    // If a column sum or a row sum is zero, the table is non-sense
152    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
153      std::cerr << "Warning: Fisher: Table is not valid\n";
154      return 1;
155    }
156
157    return oddsratio(a_,b_,c_,d_); 
158  }
159 
160  double Fisher::score(const classifier::Target& target, 
161                       const classifier::DataLookupWeighted1D& value) 
162  {
163    weighted_=true;
164    a_=b_=c_=d_=0;
165    for (size_t i=0; i<target.size(); i++)
166      if (target.binary(i))
167        if (value.data(i)>value_cutoff_)
168          a_+=value.weight(i);
169        else
170          c_+=value.weight(i);
171      else
172        if (value.data(i)>value_cutoff_)
173          b_+=value.weight(i);
174        else
175          d_+=value.weight(i);
176       
177    // If a column sum or a row sum is zero, the table is non-sense
178    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
179      // Peter should throw an exception here
180      std::cerr << "Warning: Fisher: Table is not valid\n";
181      return 1;
182    }
183
184    return oddsratio(a_,b_,c_,d_); 
185  }
186 
187    double Fisher::score(const classifier::Target& target, 
188                         const utility::vector& value, 
189                         const utility::vector& weight) 
190  {
191    weighted_=true;
192    a_=b_=c_=d_=0;
193    for (size_t i=0; i<target.size(); i++)
194      if (target.binary(i))
195        if (value(i)>value_cutoff_)
196          a_+=weight(i);
197        else
198          c_+=weight(i);
199      else
200        if (value(i)>value_cutoff_)
201          b_+=weight(i);
202        else
203          d_+=weight(i);
204       
205    // If a column sum or a row sum is zero, the table is non-sense
206    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
207      // Peter should throw an exception
208      std::cerr << "Warning: Fisher: Table is not valid\n";
209      return 1;
210    }
211
212    return oddsratio(a_,b_,c_,d_); 
213  }
214 
215  double Fisher::score(const u_int a, const u_int b, 
216                       const u_int c, const u_int d) 
217  {
218    a_=a;
219    b_=b;
220    c_=c;
221    d_=d;
222    return oddsratio(a,b,c,d); 
223  }
224
225}} // of namespace statistics and namespace theplu
Note: See TracBrowser for help on using the repository browser.