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

// $Id: Fisher.cc 675 2006-10-10 12:08:45Z jari $

/*
  Copyright (C) The authors contributing to this file.

  This file is part of the yat library, http://lev.thep.lu.se/trac/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "yat/statistics/Fisher.h"
#include "yat/statistics/Score.h"
#include "yat/statistics/utility.h"
#include "yat/classifier/DataLookupWeighted1D.h"
#include "yat/classifier/Target.h"

#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>

namespace theplu {
namespace statistics {

  
  Fisher::Fisher(bool b)
    : Score(b), a_(0), b_(0), c_(0), d_(0), oddsratio_(1.0), value_cutoff_(0)
  {
  }


  double Fisher::Chi2() const
  {
    double a,b,c,d;
    expected(a,b,c,d);
    return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b + 
      (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
  }


  void Fisher::expected(double& a, double& b, double& c, double& d) const
  {
    double N = a_+b_+c_+d_;
    a =((a_+b_)*(a_+c_)) / N;
    b =((a_+b_)*(b_+d_)) / N;
    c =((c_+d_)*(a_+c_)) / N;
    d =((c_+d_)*(b_+d_)) / N;
  }

  double Fisher::oddsratio(const double a,
                           const double b,
                           const double c,
                           const double d) 
  {
    // If a column sum or a row sum is zero, the table is nonsense
    if ((a==0 || d==0) && (c==0 || b==0)){
      //Peter, should throw exception
      std::cerr << "Warning: Fisher: Table is not valid\n";
      return oddsratio_ = 1.0;
    }

    oddsratio_=(a*d)/(b*d);
    if (absolute_ && oddsratio_<1)
      return 1/oddsratio_;
    return oddsratio_;
  }


  double Fisher::p_value() const
  {
    double p=1;
    if ( ( a_<minimum_size_ || b_<minimum_size_ || 
           c_<minimum_size_ || d_<minimum_size_) && !weighted_ ){

      p=p_value_exact();
    }
    else
      p=p_value_approximative();

    if (!absolute_){
      p=p/2;
      if (oddsratio_<0.5){
        // last term because >= not equal to !(<=)
        u_int a = static_cast<u_int>(a_);
        u_int b = static_cast<u_int>(b_);
        u_int c = static_cast<u_int>(c_);
        u_int d = static_cast<u_int>(d_);
        return 1-p+gsl_ran_hypergeometric_pdf(a, a+b, c+d, a+c);
      }
    }
    return p;
  }


  double Fisher::p_value_approximative() const
  {
    return gsl_cdf_chisq_Q(Chi2(), 1.0);
  }

  double Fisher::p_value_exact() const
  { 
    u_int a = static_cast<u_int>(a_);
    u_int b = static_cast<u_int>(b_);
    u_int c = static_cast<u_int>(c_);
    u_int d = static_cast<u_int>(d_);
    
    // Since the calculation is symmetric and cdf_hypergeometric_P
    // loops up to k we choose the smallest number to be k and mirror
    // the matrix. This choice makes the p-value two-sided.

    if (a<b && a<c && a<d)
      return statistics::cdf_hypergeometric_P(a,a+b,c+d,a+c);
    else if (b<a && b<c && b<d)
      return  statistics::cdf_hypergeometric_P(b,a+b,c+d,b+d);
    else if (c<a && c<b && c<d)
      return  statistics::cdf_hypergeometric_P(c,c+d,a+b,a+c);

    return statistics::cdf_hypergeometric_P(d,c+d,a+b,b+d);

  }


    double Fisher::score(const classifier::Target& target, 
                         const utility::vector& value)
  {
    weighted_=false;
    a_=b_=c_=d_=0;
    for (size_t i=0; i<target.size(); i++)
      if (target.binary(i))
        if (value(i)>value_cutoff_)
          a_++;
        else
          c_++;
      else
        if (value(i)>value_cutoff_)
          b_++;
        else
          d_++;
        
    // If a column sum or a row sum is zero, the table is non-sense
    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
      std::cerr << "Warning: Fisher: Table is not valid\n";
      return 1;
    }

    return oddsratio(a_,b_,c_,d_); 
  }
  
  double Fisher::score(const classifier::Target& target, 
                       const classifier::DataLookupWeighted1D& value) 
  {
    weighted_=true;
    a_=b_=c_=d_=0;
    for (size_t i=0; i<target.size(); i++)
      if (target.binary(i))
        if (value.data(i)>value_cutoff_)
          a_+=value.weight(i);
        else
          c_+=value.weight(i);
      else
        if (value.data(i)>value_cutoff_)
          b_+=value.weight(i);
        else
          d_+=value.weight(i);
        
    // If a column sum or a row sum is zero, the table is non-sense
    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
      // Peter should throw an exception here
      std::cerr << "Warning: Fisher: Table is not valid\n";
      return 1;
    }

    return oddsratio(a_,b_,c_,d_); 
  }
  
    double Fisher::score(const classifier::Target& target, 
                         const utility::vector& value, 
                         const utility::vector& weight) 
  {
    weighted_=true;
    a_=b_=c_=d_=0;
    for (size_t i=0; i<target.size(); i++)
      if (target.binary(i))
        if (value(i)>value_cutoff_)
          a_+=weight(i);
        else
          c_+=weight(i);
      else
        if (value(i)>value_cutoff_)
          b_+=weight(i);
        else
          d_+=weight(i);
        
    // If a column sum or a row sum is zero, the table is non-sense
    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
      // Peter should throw an exception
      std::cerr << "Warning: Fisher: Table is not valid\n";
      return 1;
    }

    return oddsratio(a_,b_,c_,d_); 
  }
  
  double Fisher::score(const u_int a, const u_int b, 
                       const u_int c, const u_int d) 
  {
    a_=a;
    b_=b;
    c_=c;
    d_=d;
    return oddsratio(a,b,c,d); 
  }

}} // of namespace statistics and namespace theplu