Changeset 2718

Timestamp:

Apr 12, 2012, 2:55:53 AM (11 years ago)

Author:

Peter

Message:

implement ROC p-value handling weights and ties. closes #211

Location:

trunk

Files:

• test/roc.cc (modified) (4 diffs)
• yat/statistics/ROC.cc (modified) (5 diffs)
• yat/statistics/ROC.h (modified) (7 diffs)

trunk/test/roc.cc

@@ -53 +53 @@
 void test_p_exact_with_ties(test::Suite& suite);
 void test_p_with_weights(test::Suite& suite);
+void test_p_with_weights_and_ties(test::Suite& suite);
 void test_ties(test::Suite& suite);
 
@@ -246 +247 @@
   test_p_exact_weighted(suite);
   test_p_approx_weighted(suite);
+  test_p_with_weights_and_ties(suite);
+}
+
+
+void test_p_with_weights_and_ties(test::Suite& suite)
+{
+  suite.out() << "test p with weights and ties\n";
+  statistics::ROC roc;
+  roc.add(10, true, 1.0);
+  roc.add(10, false, 3.0);
+  roc.add(20, true, 2.0);
+  roc.add(30, true, 1.0);
+  if (!suite.equal(roc.area(), 0.875)) {
+    suite.add(false);
+    suite.out() << "roc area: " << roc.area() << "\n";
+  }
+  double p = roc.p_value_one_sided();
+  if (!suite.equal(p, 8.0/24.0)) {
+    suite.add(false);
+    suite.out() << "p_value_one_sided() failed\n";
+  }
+  p = roc.p_value();
+  if (!suite.equal(p, (8.0+6.0)/24.0)) {
+    suite.add(false);
+    suite.out() << "p_value() failed\n";
+  }
 }
 
@@ -306 +333 @@
   unsigned long perm = 0;
   unsigned long k = 0;
+  unsigned long k2 = 0;
   while (true) {
     ++perm;
@@ -313 +341 @@
     if (roc2.area() >= roc.area())
       ++k;
+    if (roc2.area() <= 1-roc.area()+1e-10)
+      ++k2;
 
     if (!next_permutation(w.begin(), w.end()))
       break;
   }
-  if (!suite.xadd(suite.equal(roc.p_value_one_sided(),
-                              static_cast<double>(k)/perm))) {
+  double p_value = roc.p_value_one_sided();
+  roc.p_value_one_sided();
+  if (!suite.add(suite.equal(p_value, static_cast<double>(k)/perm))) {
     suite.out() << "area: " << roc.area() << "\n"
                 << perm << " permutations of which\n"
                 << k << " with larger (or equal) area "
                 << "corresponding to P=" << static_cast<double>(k)/perm << "\n"
-                << "p_value_one_sided() returned: " << roc.p_value_one_sided()
+                << "p_value_one_sided() returned: " << p_value
                 << "\n";
   }
-
+  p_value = roc.p_value();
+  if (!suite.add(suite.equal(p_value, static_cast<double>(k+k2)/perm))) {
+    suite.out() << "area: " << roc.area() << "\n"
+                << perm << " permutations of which\n"
+                << k << " with larger (or equal) area and\n"
+                << k2 << " with smaller (or equal) area\n"
+                << "corresponding to P="
+                << static_cast<double>(k+k2)/perm << "\n"
+                << "p_value() returned: " << p_value
+                << "\n";
+  }
 }
 

trunk/yat/statistics/ROC.cc

@@ -110 +110 @@
 
 
+  bool ROC::is_weighted(void) const
+  {
+    return pos_weights_.variance() || neg_weights_.variance()
+      || pos_weights_.mean() != neg_weights_.mean();
+  }
+
   unsigned int& ROC::minimum_size(void)
   {
@@ -146 +152 @@
 
 
-  double ROC::p_exact(double area) const
-  {
+  double ROC::p_exact_left(double area) const
+  {
+    if (is_weighted())
+      return p_left_weighted(area);
+    return p_exact_with_ties(multimap_.rbegin(), multimap_.rend(),
+                             (1-area)*pos_weights_.n()*neg_weights_.n(),
+                             pos_weights_.n(), neg_weights_.n());
+  }
+
+
+  double ROC::p_exact_right(double area) const
+  {
+    if (is_weighted())
+      return p_right_weighted(area);
     return p_exact_with_ties(multimap_.begin(), multimap_.end(),
                              area*pos_weights_.n()*neg_weights_.n(),
                              pos_weights_.n(), neg_weights_.n());
+  }
+
+
+  double ROC::p_left_weighted(double area) const
+  {
+    return count(utility::pair_first_iterator(multimap_.begin()),
+                 utility::pair_first_iterator(multimap_.end()), 1-area);
+  }
+
+
+  double ROC::p_right_weighted(double area) const
+  {
+    return count(utility::pair_first_iterator(multimap_.rbegin()),
+                 utility::pair_first_iterator(multimap_.rend()), area);
   }
 
@@ -166 +198 @@
       double p = 0;
       double abs_area = std::max(area, 1-area);
-      p = p_exact(abs_area);
+      p = p_exact_right(abs_area);
       if (has_ties_) {
-        p += p_exact_with_ties(multimap_.rbegin(), multimap_.rend(),
-                               abs_area*pos_weights_.n()*neg_weights_.n(),
-                               pos_weights_.n(), neg_weights_.n());
+        p += p_exact_left(1.0 - abs_area);
       }
       else
@@ -191 +221 @@
       return std::numeric_limits<double>::quiet_NaN();
     if (use_exact_method())
-      return p_exact(area);
+      return p_exact_right(area);
     return get_p_approx(area);
   }
@@ -211 +241 @@
   }
 
+
+  ROC::Weights::Weights(void)
+    : small_pos(0), small_neg(0), tied_pos(0), tied_neg(0)
+  {}
+
 }}} // of namespace statistics, yat, and theplu

trunk/yat/statistics/ROC.h

@@ -7 +7 @@
   Copyright (C) 2004 Peter Johansson
   Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
-  Copyright (C) 2011 Peter Johansson
+  Copyright (C) 2011, 2012 Peter Johansson
 
   This file is part of the yat library, http://dev.thep.lu.se/yat
@@ -26 +26 @@
 
 #include "Averager.h"
+#include "yat/utility/stl_utility.h"
 #include "yat/utility/yat_assert.h"
 
@@ -133 +134 @@
        \b Exact \b method: In the exact method the function goes
        through all permutations and counts what fraction for which the
-       area is greater (or equal) than area in original permutation.
+       area is greater (or equal) than area in original
+       permutation. In case all non-zero weights are not equal,
+       iterating through all permutations is not sufficient so
+       algorithm goes through all combinations instead which quickly
+       becomes a large number (N!).
 
        \b Large-sample \b Approximation: When many data points are
@@ -194 +199 @@
     typedef std::multimap<double, std::pair<bool, double> > Map;
 
+    // struct used i count functions
+    struct Weights
+    {
+      Weights(void);
+      double small_pos;
+      double small_neg;
+      double tied_pos;
+      double tied_neg;
+    };
+
     /// Implemented as in MatLab 13.1
     double get_p_approx(double) const;
 
     /**
+       return false if all non-zero weights are equal
+     */
+    bool is_weighted(void) const;
+
+    /**
        return (sum x)^2 / sum x^2
      */
@@ -203 +223 @@
 
     /*
+      Calculate probability to get an area equal (smaller) than \a
+      area given the distribution of weights and ties in multimap_
+     */
+    double p_left_weighted(double area) const;
+
+    /*
+      Calculate probability to get an area equal (greater) than \a
+      area given the distribution of weights and ties in multimap_
+     */
+    double p_right_weighted(double area) const;
+
+    /*
+      Count number of combinations (of N!) that gives weight sum equal
+      or larger than \a threshold.
+
+      Range [first, last) is used to check for ties. If, e.g., *first
+      and *(first+1) are equal implies that the two largest values are
+      equal.
+     */
+    template <typename Iterator>
+    double count(Iterator first, Iterator last, double threshold) const;
+
+    /*
+      Loop over all elements in \a weights and call count(7)
+     */
+    template <typename Iterator>
+    double count(Map& weights, Iterator iter, Iterator last,
+                 double threshold, double sum, const Weights& weight) const;
+
+    /*
+      Count number of combinations in which sum>=threshold given
+      classes and weights in \a weight. Range [iter, last) is used to
+      handle ties.
+     */
+    template <typename Iterator>
+    double count(Map& weights, Iterator iter, Iterator last,
+                 double threshold, double sum, Weights weight,
+                 const std::pair<bool, double>& entry) const;
+
+    /*
+      Calculates probability to get \a block number of pairs correctly
+      sorted when having \a pos positive samples and \a neg negative
+      samples given the distribution of ties as in [first, last).
      */
     template<typename ForwardIterator>
@@ -210 +273 @@
 
     /**
-       \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;
+       \return P(auc >= area)
+     */
+    double p_exact_right(double area) const;
+
+    /**
+       \return P(auc <= area)
+     */
+    double p_exact_left(double area) const;
 
     bool use_exact_method(void) const;
@@ -273 +337 @@
   }
 
+
+  template <typename Iterator>
+  double ROC::count(Iterator first, Iterator last, double threshold) const
+  {
+    Map map(multimap_);
+    ROC::Weights w;
+    w.small_pos = pos_weights_.sum_x();
+    w.small_neg = neg_weights_.sum_x();
+    return count(map, first, last, threshold*w.small_pos*w.small_neg, 0, w);
+  }
+
+
+
+  template <typename Iterator>
+  double ROC::count(Map& weights, Iterator iter, Iterator last,
+                    double threshold, double sum, const Weights& w) const
+  {
+    double result = 0.0;
+    // loop over all elements
+    for (Map::iterator i=weights.begin(); i!=weights.end(); ++i) {
+      Map::value_type save = *i;
+      Map::iterator hint = i;
+      ++hint;
+      weights.erase(i);
+      result += count(weights, iter, last, threshold, sum, w, save.second);
+      i = weights.insert(hint, save);
+    }
+    YAT_ASSERT(weights.size());
+    return result/weights.size();
+  }
+
+  template <typename Iterator>
+  double ROC::count(Map& weights, Iterator iter, Iterator last,
+                    double threshold, double sum, Weights w,
+                    const std::pair<bool, double>& entry) const
+  {
+    double tiny = 10e-10;
+
+    Iterator next(iter);
+    ++next;
+
+    // update weights
+    if (entry.first) {
+      w.tied_pos += entry.second;
+      w.small_pos -= entry.second;
+    }
+    else {
+      w.tied_neg += entry.second;
+      w.small_neg -= entry.second;
+    }
+
+    // last entry in equal range
+    if (next==last || *next!=*iter) {
+      sum += 0.5*w.tied_pos*w.tied_neg + w.tied_pos * w.small_neg;
+      w.tied_pos=0;
+      w.tied_neg=0;
+    }
+
+    // max sum happens if all pos values belong to current equal range
+    // and none of the remaining neg values
+    double max_sum = sum + 0.5*(w.tied_pos+w.small_pos)*w.tied_neg +
+      (w.tied_pos+w.small_pos)*w.small_neg;
+
+    if (max_sum<threshold-tiny)
+      return 0.0;
+    if (sum >= threshold-tiny)
+      return 1.0;
+
+    if (next!=last)
+      return count(weights, next, last, threshold, sum, w);
+    return 0.0;
+  }
+
 }}} // of namespace statistics, yat, and theplu
 #endif