source: trunk/lib/statistics/Fisher.cc @ 295

// $Id: Fisher.cc 295 2005-04-29 09:15:58Z peter $

#include <c++_tools/statistics/Fisher.h>
#include <c++_tools/statistics/Score.h>
#include <c++_tools/statistics/Statistics.h>

namespace theplu {
namespace statistics {

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


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

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


  double Fisher::p_value() const
  {
    // 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.
    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_);
    else 
      return statistics::cdf_hypergeometric_P(d_,c_+d_,a_+b_,b_+d_);
  }


  double Fisher::score(const gslapi::vector& x, 
                       const gslapi::vector& y, 
                       const std::vector<size_t>& train_set)
  {
    if (!train_set.size()){
      train_set_.resize(0);
      for (size_t i=0; i<target_.size(); i++)
        train_set_.push_back(i);  
    }
    else
      train_set_ = train_set;
    a_=b_=c_=d_=0;
    for (size_t i=0; i<train_set_.size(); i++)
      if (x(train_set_[i])==1)
        if (y(train_set_[i])==1)
          a_++;
        else
          c_++;
      else
        if (y(train_set_[i])==1)
          b_++;
        else
          d_++;
        
    // If a column sum or a row sum is zero, the table is nonsense
    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
      std::cerr << "Warning: Fisher: Table is not valid\n";
      return 0;
    }

    return oddsratio(a_,b_,c_,d_); 
  }
  
  double Fisher::score(const gslapi::vector& x, 
                       const gslapi::vector& y,              
                       const gslapi::vector& w, 
                       const std::vector<size_t>& train_set)
  {
    if (!train_set.size()){
      train_set_.resize(0);
      for (size_t i=0; i<target_.size(); i++)
        train_set_.push_back(i);  
    }
    else
      train_set_ = train_set;

    // Since table elements might be non-integer, we choose to not
    // store the table as member variables, which means p_value
    // function can not be used (which makes sense since p_value for
    // weighted version is not defined).

    double a=0;
    double b=0;
    double c=0;
    double d=0;
    for (size_t i=0; i<train_set_.size(); i++)
      if (x(train_set_[i])==1)
        if (y(train_set_[i])==1)
          a+=w(train_set_[i]);
        else
          c+=w(train_set_[i]);
      else
        if (y(train_set_[i])==1)
          b+=w(train_set_[i]);
        else
          d+=w(train_set_[i]);

    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