Changeset 2708

Timestamp:

Mar 14, 2012, 8:06:01 AM (11 years ago)

Author:

Peter

Message:

tests for weighted ROC p-value. refs #211

File:

• trunk/test/roc.cc (modified) (4 diffs)

trunk/test/roc.cc

@@ -24 +24 @@
 
 #include "yat/classifier/DataLookupWeighted1D.h"
+#include "yat/random/random.h"
 #include "yat/classifier/Target.h"
+#include "yat/statistics/Averager.h"
 #include "yat/statistics/Fisher.h"
 #include "yat/statistics/ROC.h"
@@ -30 +32 @@
 #include "yat/utility/Vector.h"
 
+#include <gsl/gsl_cdf.h>
+
 #include <cassert>
 #include <cmath>
+#include <deque>
 #include <fstream>
 #include <iostream>
-
+#include <set>
+#include <vector>
 
 using namespace theplu::yat;
 
 void test_empty(test::Suite&);
+void test_p_approx_weighted(test::Suite& suite);
+void test_p_approx_with_ties(test::Suite& suite);
+void test_p_approx(test::Suite& suite);
+void test_p_double_weights(test::Suite& suite);
+void test_p_exact(test::Suite& suite);
+void test_p_exact_weighted(test::Suite& suite);
+void test_p_exact_with_ties(test::Suite& suite);
+void test_p_with_weights(test::Suite& suite);
 void test_ties(test::Suite& suite);
-void test_p_exact(test::Suite& suite);
-void test_p_approx(test::Suite& suite);
-void test_p_exact_with_ties(test::Suite& suite);
-void test_p_approx_with_ties(test::Suite& suite);
 
 int main(int argc, char* argv[])
@@ -108 +118 @@
   test_p_exact(suite);
   test_empty(suite);
+  test_p_with_weights(suite);
   return suite.return_value();
 }
@@ -227 +238 @@
   suite.err() << "test empty done\n";
 }
+
+
+void test_p_with_weights(test::Suite& suite)
+{
+  suite.out() << "test p with weights\n";
+  test_p_double_weights(suite);
+  test_p_exact_weighted(suite);
+  test_p_approx_weighted(suite);
+}
+
+
+void test_p_double_weights(test::Suite& suite)
+{
+  suite.out() << "test p with double weights\n";
+  std::vector<double> w(5,1.0);
+  w[0]=0.1;
+  w[4]=10;
+  std::vector<double> x(5);
+  for (size_t i=0; i<x.size(); ++i)
+    x[i] = i;
+  statistics::ROC roc;
+  statistics::ROC roc2;
+  for (size_t i=0; i<x.size(); ++i) {
+    roc.add(x[i], i<2, w[i]);
+    roc2.add(x[i], i<2, 2*w[i]);
+  }
+  if (!suite.equal(roc.area(), roc2.area())) {
+    suite.add(false);
+    suite.err() << "area failed\n";
+  }
+  if (!suite.equal(roc.p_value_one_sided(), roc2.p_value_one_sided())) {
+    suite.add(false);
+    suite.err() << "exact p failed\n";
+  }
+  roc.minimum_size() = 0;
+  roc2.minimum_size() = 0;
+  if (!suite.equal(roc.p_value_one_sided(), roc2.p_value_one_sided())) {
+    suite.xadd(false);
+    suite.err() << "approximative p failed\n";
+  }
+}
+
+
+void test_p_exact_weighted(test::Suite& suite)
+{
+  suite.out() << "test p exact weighted\n";
+  std::vector<double> x(10);
+  std::vector<std::pair<bool, double> > w(10);
+  for (size_t i=0; i<x.size(); ++i) {
+    x[i] = i;
+    w[i].first = false;
+    w[i].second = 1.0;
+  }
+  w[3].first = true;
+  w[7].first = true;
+  w[8].first = true;
+  w[1].second = 10;
+  w[7].second = 10;
+  w[9].second = 0.1;
+
+  statistics::ROC roc;
+  for (size_t i=0; i<x.size(); ++i)
+    roc.add(x[i], w[i].first, w[i].second);
+  roc.minimum_size() = 0;
+
+  std::sort(w.begin(), w.end());
+  unsigned long perm = 0;
+  unsigned long k = 0;
+  while (true) {
+    ++perm;
+    statistics::ROC roc2;
+    for (size_t i=0; i<x.size(); ++i)
+      roc2.add(x[i], w[i].first, w[i].second);
+    if (roc2.area() >= roc.area())
+      ++k;
+
+    if (!next_permutation(w.begin(), w.end()))
+      break;
+  }
+  if (!suite.xadd(suite.equal(roc.p_value_one_sided(),
+                              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()
+                << "\n";
+  }
+
+}
+
+
+void test_p_approx_weighted(test::Suite& suite)
+{
+  suite.out() << "test p approx weighted\n";
+  std::vector<double> x(100);
+  std::vector<double> w(100, 1.0);
+  std::deque<bool> label(100);
+
+  for (size_t i=0; i<x.size(); ++i) {
+    x[i] = i;
+    label[i] = i>30 && i<70;
+    w[i] = 100.0 / (i+100);
+  }
+
+  statistics::ROC roc;
+  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;
+  checkpoints.insert(10);
+  checkpoints.insert(100);
+  checkpoints.insert(1000);
+  checkpoints.insert(10000);
+  checkpoints.insert(100000);
+  checkpoints.insert(perm);
+  statistics::Averager averager;
+  for (size_t i=1; i<=perm; ++i) {
+    random::random_shuffle(x.begin(), x.end());
+    statistics::ROC roc2;
+    for (size_t j=0; j<x.size(); ++j)
+      roc2.add(x[j], label[j], w[j]);
+    if (roc2.area()>=roc.area())
+      averager.add(1.0);
+    else
+      averager.add(0.0);
+    if (checkpoints.find(i)!=checkpoints.end()) {
+      if (gsl_cdf_binomial_P(averager.sum_x(), p, averager.n())<1e-10 ||
+          gsl_cdf_binomial_Q(averager.sum_x(), p, averager.n())<1e-10) {
+        suite.err() << "error: approx p value and permutation p-value "
+                    << "deviate more than expected\n"
+                    << "approx p: " << p << "\n"
+                    << "permutations: " << averager.n() << "\n"
+                    << "successful: " << averager.sum_x() << "\n"
+                    << "corresponds to P=" << averager.mean() << "\n";
+        suite.xadd(false);
+        return;
+      }
+    }
+  }
+}