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

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

Addresses #153. Introduced yat namespace. Removed alignment namespace. Clean up of code.

  • 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 680 2006-10-11 17:49:03Z 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  double Fisher::oddsratio(const double a,
63                           const double b,
64                           const double c,
65                           const double d) 
66  {
67    // If a column sum or a row sum is zero, the table is nonsense
68    if ((a==0 || d==0) && (c==0 || b==0)){
69      //Peter, should throw exception
70      std::cerr << "Warning: Fisher: Table is not valid\n";
71      return oddsratio_ = 1.0;
72    }
73
74    oddsratio_=(a*d)/(b*d);
75    if (absolute_ && oddsratio_<1)
76      return 1/oddsratio_;
77    return oddsratio_;
78  }
79
80
81  double Fisher::p_value() const
82  {
83    double p=1;
84    if ( ( a_<minimum_size_ || b_<minimum_size_ || 
85           c_<minimum_size_ || d_<minimum_size_) && !weighted_ ){
86
87      p=p_value_exact();
88    }
89    else
90      p=p_value_approximative();
91
92    if (!absolute_){
93      p=p/2;
94      if (oddsratio_<0.5){
95        // last term because >= not equal to !(<=)
96        u_int a = static_cast<u_int>(a_);
97        u_int b = static_cast<u_int>(b_);
98        u_int c = static_cast<u_int>(c_);
99        u_int d = static_cast<u_int>(d_);
100        return 1-p+gsl_ran_hypergeometric_pdf(a, a+b, c+d, a+c);
101      }
102    }
103    return p;
104  }
105
106
107  double Fisher::p_value_approximative() const
108  {
109    return gsl_cdf_chisq_Q(Chi2(), 1.0);
110  }
111
112  double Fisher::p_value_exact() const
113  { 
114    u_int a = static_cast<u_int>(a_);
115    u_int b = static_cast<u_int>(b_);
116    u_int c = static_cast<u_int>(c_);
117    u_int d = static_cast<u_int>(d_);
118   
119    // Since the calculation is symmetric and cdf_hypergeometric_P
120    // loops up to k we choose the smallest number to be k and mirror
121    // the matrix. This choice makes the p-value two-sided.
122
123    if (a<b && a<c && a<d)
124      return statistics::cdf_hypergeometric_P(a,a+b,c+d,a+c);
125    else if (b<a && b<c && b<d)
126      return  statistics::cdf_hypergeometric_P(b,a+b,c+d,b+d);
127    else if (c<a && c<b && c<d)
128      return  statistics::cdf_hypergeometric_P(c,c+d,a+b,a+c);
129
130    return statistics::cdf_hypergeometric_P(d,c+d,a+b,b+d);
131
132  }
133
134
135    double Fisher::score(const classifier::Target& target, 
136                         const utility::vector& value)
137  {
138    weighted_=false;
139    a_=b_=c_=d_=0;
140    for (size_t i=0; i<target.size(); i++)
141      if (target.binary(i))
142        if (value(i)>value_cutoff_)
143          a_++;
144        else
145          c_++;
146      else
147        if (value(i)>value_cutoff_)
148          b_++;
149        else
150          d_++;
151       
152    // If a column sum or a row sum is zero, the table is non-sense
153    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
154      std::cerr << "Warning: Fisher: Table is not valid\n";
155      return 1;
156    }
157
158    return oddsratio(a_,b_,c_,d_); 
159  }
160 
161  double Fisher::score(const classifier::Target& target, 
162                       const classifier::DataLookupWeighted1D& value) 
163  {
164    weighted_=true;
165    a_=b_=c_=d_=0;
166    for (size_t i=0; i<target.size(); i++)
167      if (target.binary(i))
168        if (value.data(i)>value_cutoff_)
169          a_+=value.weight(i);
170        else
171          c_+=value.weight(i);
172      else
173        if (value.data(i)>value_cutoff_)
174          b_+=value.weight(i);
175        else
176          d_+=value.weight(i);
177       
178    // If a column sum or a row sum is zero, the table is non-sense
179    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
180      // Peter should throw an exception here
181      std::cerr << "Warning: Fisher: Table is not valid\n";
182      return 1;
183    }
184
185    return oddsratio(a_,b_,c_,d_); 
186  }
187 
188    double Fisher::score(const classifier::Target& target, 
189                         const utility::vector& value, 
190                         const utility::vector& weight) 
191  {
192    weighted_=true;
193    a_=b_=c_=d_=0;
194    for (size_t i=0; i<target.size(); i++)
195      if (target.binary(i))
196        if (value(i)>value_cutoff_)
197          a_+=weight(i);
198        else
199          c_+=weight(i);
200      else
201        if (value(i)>value_cutoff_)
202          b_+=weight(i);
203        else
204          d_+=weight(i);
205       
206    // If a column sum or a row sum is zero, the table is non-sense
207    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
208      // Peter should throw an exception
209      std::cerr << "Warning: Fisher: Table is not valid\n";
210      return 1;
211    }
212
213    return oddsratio(a_,b_,c_,d_); 
214  }
215 
216  double Fisher::score(const u_int a, const u_int b, 
217                       const u_int c, const u_int d) 
218  {
219    a_=a;
220    b_=b;
221    c_=c;
222    d_=d;
223    return oddsratio(a,b,c,d); 
224  }
225
226}}} // of namespace statistics, yat and theplu
Note: See TracBrowser for help on using the repository browser.