Ignore:
Timestamp:
Mar 24, 2013, 1:51:14 AM (10 years ago)
Author:
Peter
Message:

refs #689. Deprecate Fisher::one_sided_p; implement Fisher::left_p and right_p.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/statistics/Fisher.cc

    r2919 r3004  
    55  Copyright (C) 2005 Peter Johansson
    66  Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
    7   Copyright (C) 2009, 2010, 2011, 2012 Peter Johansson
     7  Copyright (C) 2009, 2010, 2011, 2012, 2013 Peter Johansson
    88
    99  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    4040namespace statistics {
    4141
    42  
     42
    4343  Fisher::Fisher()
    4444    : a_(0), b_(0), c_(0), d_(0), minimum_size_(10), oddsratio_(1.0)
     
    5454  bool Fisher::calculate_p_exact(void) const
    5555  {
    56     return ( a_<minimum_size_ || b_<minimum_size_ || 
    57              c_<minimum_size_ || d_<minimum_size_); 
    58   }
    59 
    60   double Fisher::Chi2() const
     56    return ( a_<minimum_size_ || b_<minimum_size_ ||
     57             c_<minimum_size_ || d_<minimum_size_);
     58  }
     59
     60  double Fisher::Chi2(void) const
    6161  {
    6262    double a,b,c,d;
    6363    expected(a,b,c,d);
    64     return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b + 
     64    return (a-a_)*(a-a_)/a + (b-b_)*(b-b_)/b +
    6565      (c-c_)*(c-c_)/c + (d-d_)*(d-d_)/d;
    6666  }
     
    9797                           const unsigned int b,
    9898                           const unsigned int c,
    99                            const unsigned int d) 
     99                           const unsigned int d)
    100100  {
    101101    // If a column sum or a row sum is zero, the table is nonsense
     
    108108    d_ = d;
    109109
    110     oddsratio_= (static_cast<double>(a) / static_cast<double>(b) * 
     110    oddsratio_= (static_cast<double>(a) / static_cast<double>(b) *
    111111                 static_cast<double>(d) / static_cast<double>(c) );
    112112    return oddsratio_;
     
    120120
    121121
    122   double Fisher::p_value() const
     122  double Fisher::p_left(void) const
     123  {
     124    if (!calculate_p_exact()) {
     125      if (oddsratio_>1.0)
     126        return 1.0-p_value_approximative()/2.0;
     127      return p_value_approximative()/2.0;
     128    }
     129    // check for overflow
     130    assert(c_ <= c_+d_ && d_ <= c_+d_);
     131    assert(a_ <= a_+b_ && b_ <= a_+b_);
     132    assert(a_ <= a_+c_ && c_ <= a_+c_);
     133
     134    return cdf_hypergeometric_P(a_, a_+b_, c_+d_, a_+c_);
     135  }
     136
     137
     138  double Fisher::p_value(void) const
    123139  {
    124140    if (calculate_p_exact())
     
    128144
    129145
    130   double Fisher::p_value_one_sided() const
     146  double Fisher::p_right(void) const
    131147  {
    132148    if (!calculate_p_exact()) {
     
    140156    assert(a_ <= a_+c_ && c_ <= a_+c_);
    141157
    142     return statistics::cdf_hypergeometric_P(c_, c_+d_, a_+b_, a_+c_);
    143   }
    144 
    145 
    146   double Fisher::p_value_approximative() const
     158    return cdf_hypergeometric_P(c_, c_+d_, a_+b_, a_+c_);
     159  }
     160
     161
     162  double Fisher::p_value_one_sided(void) const
     163  {
     164    return p_right();
     165  }
     166
     167
     168  double Fisher::p_value_approximative(void) const
    147169  {
    148170    return gsl_cdf_chisq_Q(Chi2(), 1.0);
    149171  }
    150172
    151   double Fisher::p_value_exact() const
    152   {
     173
     174  double Fisher::p_left_exact(void) const
     175  {
     176    return 0;
     177  }
     178
     179
     180  double Fisher::p_right_exact(void) const
     181  {
     182    return 0;
     183  }
     184
     185
     186  double Fisher::p_value_exact(void) const
     187  {
    153188    double res=0;
    154189    double a, c, tmp;
    155190    expected(a, tmp, c, tmp);
    156     // sum P(x) over x for which abs(x-E(a))>=abs(a-E(a)) 
     191    // sum P(x) over x for which abs(x-E(a))>=abs(a-E(a))
    157192
    158193    // counting left tail
Note: See TracChangeset for help on using the changeset viewer.