Changeset 2585


Ignore:
Timestamp:
Oct 25, 2011, 10:43:28 PM (10 years ago)
Author:
Peter
Message:

handle ties in exact calculation of ROC p-value. fixes #144

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/roc.cc

    r2556 r2585  
    130130  roc.area();
    131131  if (!suite.equal(roc.p_value_one_sided(), 3.0/10.0)) {
    132     suite.xadd(false);
    133     suite.out() << "p_value_one_sided: expected 0.3\n";
     132    suite.add(false);
     133    suite.out() << "  p_value_one_sided: expected 0.3\n";
    134134  }
    135135  else
    136     suite.xadd(true);
     136    suite.add(true);
    137137  if (!suite.equal(roc.p_value(), 5.0/10.0)) {
    138     suite.xadd(false);
    139     suite.out() << "p_value: expected 0.5\n";
     138    suite.add(false);
     139    suite.out() << "  (two-sided) p_value: expected 0.5\n";
    140140  }
    141141  else
    142     suite.xadd(true);
     142    suite.add(true);
    143143}
    144144
  • trunk/yat/statistics/ROC.cc

    r2572 r2585  
    4747    if (!w)
    4848      return;
    49     typedef std::multimap<double, std::pair<bool, double> > Map;
    50     Map::value_type element(x, std::make_pair(target, w));
    51     Map::iterator lower = multimap_.lower_bound(x);
     49    ROC::Map::value_type element(x, std::make_pair(target, w));
     50    ROC::Map::iterator lower = multimap_.lower_bound(x);
    5251    if (lower!=multimap_.end() && lower->first == x)
    5352      has_ties_ = true;
    5453    multimap_.insert(lower, element);
    55     //    multimap_.insert(std::make_pair(x, std::make_pair(target, w)));
    5654    if (target)
    5755      w_pos_+=w;
     
    107105  }
    108106
    109   double ROC::get_p_exact(const double block, const double nof_pos,
    110                           const double nof_neg) const
    111   {
    112     if (block <= 0.0)
    113       return 1.0;
    114     if (block > nof_neg*nof_pos)
    115       return 0.0;
    116     double p1 = get_p_exact(block-nof_neg, nof_pos-1, nof_neg);
    117     double p2 = get_p_exact(block, nof_pos, nof_neg-1);
    118     return nof_pos/(nof_pos+nof_neg)*p1 + nof_neg/(nof_pos+nof_neg)*p2;
    119   }
    120107
    121108  unsigned int& ROC::minimum_size(void)
     
    149136
    150137
     138  double ROC::p_exact(double area) const
     139  {
     140    if (area>=0.5)
     141      return p_exact_with_ties(multimap_.begin(), multimap_.end(),
     142                               area*w_pos_*w_neg_, w_pos_, w_neg_);
     143    return p_exact_with_ties(multimap_.rbegin(), multimap_.rend(),
     144                             (1.0-area)*w_pos_*w_neg_, w_pos_, w_neg_);
     145   
     146   
     147  }
     148
     149
    151150  double ROC::p_value() const
    152   {
    153     double p = 2*p_value_one_sided();
    154     return std::min(p, 2-p);
    155   }
    156 
    157 
    158   double ROC::p_value_one_sided() const
    159151  {
    160152    double area(area_);
     
    163155      area = auc.score(multimap_);
    164156    }
    165     if (n_pos() < minimum_size_ && n_neg() < minimum_size_) {
    166       // for small areas we calculate probabilitu to get larger area -
    167       // not larger or equal.
    168       if (area<0.5)
    169         return 1-get_p_exact((1-area)*n_pos()*n_neg(), n_pos(), n_neg());
    170       return get_p_exact(area*n_pos()*n_neg(), n_pos(), n_neg());
    171     }
     157    if (use_exact_method()) {
     158      double p = 0;
     159      p = p_exact(area);
     160      if (has_ties_) {
     161        p += p_exact(1-area);
     162      }
     163      else
     164        p *= 2.0;
     165      // avoid double counting when area is 0.5
     166      return std::min(p, 1.0);
     167    }
     168    return 2*get_p_approx(area);
     169  }
     170
     171
     172  double ROC::p_value_one_sided() const
     173  {
     174    double area(area_);
     175    if (area_==-1){
     176      AUC auc(false);
     177      area = auc.score(multimap_);
     178    }
     179    if (use_exact_method())
     180      return p_exact(area);
    172181    return get_p_approx(area);
    173182  }
     
    183192  }
    184193
     194
     195  bool ROC::use_exact_method(void) const
     196  {
     197    return (n_pos() < minimum_size_) && (n_neg() < minimum_size_);
     198  }
     199
    185200}}} // of namespace statistics, yat, and theplu
  • trunk/yat/statistics/ROC.h

    r2556 r2585  
    2727#include <map>
    2828#include <utility>
     29
     30#include <gsl/gsl_randist.h>
     31
     32// debug
     33#include <iostream>
    2934
    3035namespace theplu {
     
    124129
    125130  private:
     131    typedef std::multimap<double, std::pair<bool, double> > Map;
    126132   
    127133    /// Implemented as in MatLab 13.1
    128134    double get_p_approx(double) const;
    129135
    130     /// Implemented as in MatLab 13.1
    131     double get_p_exact(const double, const double, const double) const;
     136    /*
     137     */
     138    template<typename ForwardIterator>
     139    double p_exact_with_ties(ForwardIterator first, ForwardIterator last,
     140                             double block, double pos, double neg) const;
     141
     142    /**
     143       \return probability to get auc >= \a area. If area<0.5
     144       probability to auc <= area is returned
     145
     146       \note assumes all non-zero weights are equal (typically unity
     147       but not necessarily
     148    */
     149    double p_exact(double area) const;
     150
     151    bool use_exact_method(void) const;
    132152
    133153    double area_;
     
    136156    double w_neg_;
    137157    double w_pos_;
    138     // <data pair<class, weight> >
    139     typedef std::multimap<double, std::pair<bool, double> > Map;
    140158    Map multimap_;
    141159  };
    142160
     161  template<typename ForwardIterator>
     162  double
     163  ROC::p_exact_with_ties(ForwardIterator begin, ForwardIterator end,
     164                         double block, double pos, double neg) const
     165  {
     166    /**
     167    std::cerr << "block: " << block << "\n";
     168    std::cerr << "pos  : " << pos << "\n";
     169    std::cerr << "neg  : " << neg << "\n";
     170    */
     171    if (block <= 0)
     172      return 1.0;
     173    if (block > pos*neg)
     174      return 0.0;
     175
     176    ForwardIterator iter(begin);
     177    size_t n=0;
     178    while (iter!=end && iter->first == begin->first) {
     179      ++iter;
     180      ++n;
     181    }
     182    double result = 0;
     183    //std::cerr << "n: " << n << "\n";
     184    for (size_t i=0; i<=n; ++i) {
     185      double pos1 = i;
     186      double neg1 = n-i;
     187      double pos2 = pos-pos1;
     188      double neg2 = neg-neg1;
     189      result += gsl_ran_hypergeometric_pdf(i, pos, neg, n)
     190        * p_exact_with_ties(iter, end,
     191                            block - pos2*neg1 - 0.5*pos1*neg1,
     192                            pos2, neg2);
     193      //std::cerr << "= " << n << " " << i << " " << result << "\n";
     194    }
     195    return result;
     196  }
     197
    143198}}} // of namespace statistics, yat, and theplu
    144199
Note: See TracChangeset for help on using the changeset viewer.