Changeset 2709

Timestamp:

Mar 15, 2012, 4:26:47 AM (11 years ago)

Author:

Peter

Message:

refs #211. Changed how number of samples, n, is calculated in large
sample approximation p-value calculation. Before the sum_w was used
and now sum(w)*sum(w)/sum(ww) is used.

Location:

trunk

Files:

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

trunk/test/roc.cc

@@ -275 +275 @@
   roc2.minimum_size() = 0;
   if (!suite.equal(roc.p_value_one_sided(), roc2.p_value_one_sided())) {
-    suite.xadd(false);
+    suite.add(false);
     suite.err() << "approximative p failed\n";
   }
@@ -333 +333 @@
 {
   suite.out() << "test p approx weighted\n";
-  std::vector<double> x(100);
-  std::vector<double> w(100, 1.0);
-  std::deque<bool> label(100);
+  std::vector<double> x(200);
+  std::vector<double> w(200, 1.0);
+  std::deque<bool> label(200);
 
   for (size_t i=0; i<x.size(); ++i) {
     x[i] = i;
     label[i] = i>30 && i<70;
-    w[i] = 100.0 / (i+100);
+    if (i<100)
+      w[i] = 100.0 / (100+i);
+    else
+      w[i] = 0.0001;
   }
 
@@ -346 +349 @@
   for (size_t i=0; i<x.size(); ++i)
     roc.add(x[i], label[i], w[i]);
-  suite.out() << "area: " << roc.area() << std::endl;
   roc.minimum_size() = 0;
   double p = roc.p_value_one_sided();
-  suite.out() << "p: " << p << std::endl;
 
   std::set<size_t> checkpoints;
-  size_t perm = 1000000;
+  size_t perm = 100000;
   checkpoints.insert(10);
   checkpoints.insert(100);
   checkpoints.insert(1000);
   checkpoints.insert(10000);
-  checkpoints.insert(100000);
   checkpoints.insert(perm);
   statistics::Averager averager;
@@ -374 +374 @@
         suite.err() << "error: approx p value and permutation p-value "
                     << "deviate more than expected\n"
+                    << "area: " << roc.area() << "\n"
                     << "approx p: " << p << "\n"
                     << "permutations: " << averager.n() << "\n"
                     << "successful: " << averager.sum_x() << "\n"
                     << "corresponds to P=" << averager.mean() << "\n";
-        suite.xadd(false);
+        suite.add(false);
         return;
       }

trunk/yat/statistics/ROC.cc

@@ -4 +4 @@
   Copyright (C) 2004, 2005 Peter Johansson
   Copyright (C) 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
@@ -54 +54 @@
     multimap_.insert(lower, element);
     if (target)
-      w_pos_+=w;
+      pos_weights_.add(w);
     else
-      w_neg_+=w;
+      neg_weights_.add(w);
     area_ = std::numeric_limits<double>::quiet_NaN();
   }
@@ -73 +73 @@
   double ROC::get_p_approx(double x) const
   {
+    size_t n_pos = nof_points(pos_weights_);
+    size_t n_neg = nof_points(neg_weights_);
+    size_t nof_samples = n_pos + n_neg;
     // make x standard normal
     x -= 0.5;
     // Not integrating from the middle of the bin, but from the inner edge.
     if (x>0)
-      x -= 0.5/(n_pos()*n_neg());
+      x -= 0.5/(n_pos*n_neg);
     else if(x<0)
-      x += 0.5/(n_pos()*n_neg());
+      x += 0.5/(n_pos*n_neg);
     else
       return 0.5;
-    double var = 1.0+n();
+    double var = 1.0+nof_samples;
     if (has_ties_) {
       double correction = 0;
@@ -100 +103 @@
         mn/12 [ N+1 - 1/(N(N-1)) * sum(t(t-1)(t+1)) ]
       */
-      var -= correction/(n() * (n()-1));
-    }
-    var = var / (12*n_pos()*n_neg());
+      var -= correction/(nof_samples * (nof_samples-1));
+    }
+    var = var / (12*n_pos*n_neg);
     return gsl_cdf_gaussian_Q(x, std::sqrt(var));
   }
@@ -127 +130 @@
   double ROC::n_neg(void) const
   {
-    return w_neg_;
+    return neg_weights_.sum_x();
   }
 
@@ -133 +136 @@
   double ROC::n_pos(void) const
   {
-    return w_pos_;
+    return pos_weights_.sum_x();
+  }
+
+
+  size_t ROC::nof_points(const Averager& a) const
+  {
+    return a.sum_x()*a.sum_x()/a.sum_xx();
   }
 
@@ -140 +149 @@
   {
     return p_exact_with_ties(multimap_.begin(), multimap_.end(),
-                             area*w_pos_*w_neg_, w_pos_, w_neg_);
+                             area*n_pos()*n_neg(), n_pos(), n_neg());
   }
 
@@ -159 +168 @@
       if (has_ties_) {
         p += p_exact_with_ties(multimap_.rbegin(), multimap_.rend(),
-                               abs_area*w_pos_*w_neg_, w_pos_, w_neg_);
+                               abs_area*n_pos()*n_neg(), n_pos(), n_neg());
       }
       else
@@ -189 +198 @@
     area_ = std::numeric_limits<double>::quiet_NaN();
     has_ties_ = false;
-    w_pos_=0;
-    w_neg_=0;
+    neg_weights_.reset();
+    pos_weights_.reset();
     multimap_.clear();
   }

trunk/yat/statistics/ROC.h

@@ -24 +24 @@
   along with yat. If not, see <http://www.gnu.org/licenses/>.
 */
+
+#include "Averager.h"
 
 #include <gsl/gsl_randist.h>
@@ -154 +156 @@
 
        where sum runs over different data values (of ties) and \f$ n_x
-       \f$ is number data points with that value. The sum i a
+       \f$ is number data points with that value. The sum is a
        correction term for ties and is zero if there are no ties.
 
+       The number of samples in a group, \f$ n^+ \f$, is calculated as
+       \f$ n = (\sum w)^2 / \sum w^2 \f$
+
        \return \f$ P(a \ge \textrm{area}) \f$
-
-       \note Weights should be -1, 0, or 1; otherwise the p-value is
-       undefined and may change in future versions.
      */
     double p_value_one_sided(void) const;
@@ -194 +196 @@
     double get_p_approx(double) const;
 
+    /**
+       return (sum x)^2 / sum x^2
+     */
+    size_t nof_points(const Averager& a) const;
+
     /*
      */
@@ -214 +221 @@
     bool has_ties_;
     unsigned int minimum_size_;
-    double w_neg_;
-    double w_pos_;
+    Averager neg_weights_;
+    Averager pos_weights_;
     Map multimap_;
   };