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

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

Addresses #170.

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