Changeset 2585

Timestamp:

Oct 25, 2011, 10:43:28 PM (12 years ago)

Author:

Peter

Message:

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

Location:

trunk

Files:

• test/roc.cc (modified) (1 diff)
• yat/statistics/ROC.cc (modified) (5 diffs)
• yat/statistics/ROC.h (modified) (3 diffs)

trunk/test/roc.cc

@@ -130 +130 @@
   roc.area();
   if (!suite.equal(roc.p_value_one_sided(), 3.0/10.0)) {
-    suite.xadd(false);
-    suite.out() << "p_value_one_sided: expected 0.3\n";
+    suite.add(false);
+    suite.out() << "  p_value_one_sided: expected 0.3\n";
   }
   else
-    suite.xadd(true);
+    suite.add(true);
   if (!suite.equal(roc.p_value(), 5.0/10.0)) {
-    suite.xadd(false);
-    suite.out() << "p_value: expected 0.5\n";
+    suite.add(false);
+    suite.out() << "  (two-sided) p_value: expected 0.5\n";
   }
   else
-    suite.xadd(true);
+    suite.add(true);
 }
 

trunk/yat/statistics/ROC.cc

@@ -47 +47 @@
     if (!w)
       return;
-    typedef std::multimap<double, std::pair<bool, double> > Map;
-    Map::value_type element(x, std::make_pair(target, w));
-    Map::iterator lower = multimap_.lower_bound(x);
+    ROC::Map::value_type element(x, std::make_pair(target, w));
+    ROC::Map::iterator lower = multimap_.lower_bound(x);
     if (lower!=multimap_.end() && lower->first == x)
       has_ties_ = true;
     multimap_.insert(lower, element);
-    //    multimap_.insert(std::make_pair(x, std::make_pair(target, w)));
     if (target)
       w_pos_+=w;
@@ -107 +105 @@
   }
 
-  double ROC::get_p_exact(const double block, const double nof_pos, 
-                          const double nof_neg) const
-  {
-    if (block <= 0.0)
-      return 1.0;
-    if (block > nof_neg*nof_pos)
-      return 0.0;
-    double p1 = get_p_exact(block-nof_neg, nof_pos-1, nof_neg);
-    double p2 = get_p_exact(block, nof_pos, nof_neg-1);
-    return nof_pos/(nof_pos+nof_neg)*p1 + nof_neg/(nof_pos+nof_neg)*p2;
-  }
 
   unsigned int& ROC::minimum_size(void)
@@ -149 +136 @@
 
 
+  double ROC::p_exact(double area) const
+  {
+    if (area>=0.5)
+      return p_exact_with_ties(multimap_.begin(), multimap_.end(),
+                               area*w_pos_*w_neg_, w_pos_, w_neg_);
+    return p_exact_with_ties(multimap_.rbegin(), multimap_.rend(),
+                             (1.0-area)*w_pos_*w_neg_, w_pos_, w_neg_);
+    
+   
+  }
+
+
   double ROC::p_value() const
-  {
-    double p = 2*p_value_one_sided();
-    return std::min(p, 2-p); 
-  }
-
-
-  double ROC::p_value_one_sided() const
   {
     double area(area_);
@@ -163 +155 @@
       area = auc.score(multimap_);
     }
-    if (n_pos() < minimum_size_ && n_neg() < minimum_size_) {
-      // for small areas we calculate probabilitu to get larger area -
-      // not larger or equal.
-      if (area<0.5)
-        return 1-get_p_exact((1-area)*n_pos()*n_neg(), n_pos(), n_neg());
-      return get_p_exact(area*n_pos()*n_neg(), n_pos(), n_neg());
-    }
+    if (use_exact_method()) {
+      double p = 0;
+      p = p_exact(area); 
+      if (has_ties_) {
+        p += p_exact(1-area);
+      }
+      else
+        p *= 2.0;
+      // avoid double counting when area is 0.5
+      return std::min(p, 1.0);
+    }
+    return 2*get_p_approx(area);
+  }
+
+
+  double ROC::p_value_one_sided() const
+  {
+    double area(area_);
+    if (area_==-1){
+      AUC auc(false);
+      area = auc.score(multimap_);
+    }
+    if (use_exact_method())
+      return p_exact(area);
     return get_p_approx(area);
   }
@@ -183 +192 @@
   }
 
+
+  bool ROC::use_exact_method(void) const
+  {
+    return (n_pos() < minimum_size_) && (n_neg() < minimum_size_);
+  }
+
 }}} // of namespace statistics, yat, and theplu

trunk/yat/statistics/ROC.h

@@ -27 +27 @@
 #include <map>
 #include <utility>
+
+#include <gsl/gsl_randist.h>
+
+// debug
+#include <iostream>
 
 namespace theplu {
@@ -124 +129 @@
 
   private:
+    typedef std::multimap<double, std::pair<bool, double> > Map;
     
     /// Implemented as in MatLab 13.1
     double get_p_approx(double) const;
 
-    /// Implemented as in MatLab 13.1
-    double get_p_exact(const double, const double, const double) const;
+    /*
+     */
+    template<typename ForwardIterator>
+    double p_exact_with_ties(ForwardIterator first, ForwardIterator last,
+                             double block, double pos, double neg) const;
+
+    /**
+       \return probability to get auc >= \a area. If area<0.5
+       probability to auc <= area is returned 
+
+       \note assumes all non-zero weights are equal (typically unity
+       but not necessarily
+    */
+    double p_exact(double area) const;
+
+    bool use_exact_method(void) const;
 
     double area_;
@@ -136 +156 @@
     double w_neg_;
     double w_pos_;
-    // <data pair<class, weight> >
-    typedef std::multimap<double, std::pair<bool, double> > Map;
     Map multimap_;
   };
 
+  template<typename ForwardIterator>
+  double 
+  ROC::p_exact_with_ties(ForwardIterator begin, ForwardIterator end,
+                         double block, double pos, double neg) const
+  {
+    /**
+    std::cerr << "block: " << block << "\n";
+    std::cerr << "pos  : " << pos << "\n";
+    std::cerr << "neg  : " << neg << "\n";
+    */
+    if (block <= 0)
+      return 1.0;
+    if (block > pos*neg)
+      return 0.0;
+
+    ForwardIterator iter(begin);
+    size_t n=0;
+    while (iter!=end && iter->first == begin->first) {
+      ++iter;
+      ++n;
+    }
+    double result = 0;
+    //std::cerr << "n: " << n << "\n";
+    for (size_t i=0; i<=n; ++i) {
+      double pos1 = i;
+      double neg1 = n-i;
+      double pos2 = pos-pos1;
+      double neg2 = neg-neg1;
+      result += gsl_ran_hypergeometric_pdf(i, pos, neg, n)
+        * p_exact_with_ties(iter, end, 
+                            block - pos2*neg1 - 0.5*pos1*neg1, 
+                            pos2, neg2);
+      //std::cerr << "= " << n << " " << i << " " << result << "\n";
+    }
+    return result;
+  }
+
 }}} // of namespace statistics, yat, and theplu