Changeset 1600

Timestamp:

Oct 27, 2008, 8:10:49 PM (14 years ago)

Author:

Peter

Message:

• Fixed bug that data disappeared because data internally was stored in a set (is now a multiset).
• Changed ordering of data, which implies that sign of signed_score is flipped.
• Added some tests for KolmogorovSmirnov?

Location:

trunk

Files:

• test/kolmogorov_smirnov_test.cc (modified) (1 diff)
• yat/statistics/KolmogorovSmirnov.cc (modified) (1 diff)
• yat/statistics/KolmogorovSmirnov.h (modified) (1 diff)

trunk/test/kolmogorov_smirnov_test.cc

@@ -22 +22 @@
 #include "Suite.h"
 
+#include "yat/statistics/Averager.h"
 #include "yat/statistics/KolmogorovSmirnov.h"
 
-#include <ostream>
+#include <cmath>
+#include <iostream>
+
+using namespace theplu::yat;
+
+void test_one_sample(test::Suite&);
+void test_p_value(test::Suite&);
 
 int main(int argc, char* argv[])
 {  
-  using namespace theplu::yat;
-
   test::Suite suite(argc, argv);
 
-  statistics::KolmogorovSmirnov ks;
-  for (size_t i=0; i<10; ++i)
-    ks.add(i,true);
-  for (size_t i=0; i<1000; ++i)
-    ks.add(100+i,false);
-  suite.add(suite.equal(ks.score(), 1.0));
-  for (size_t i=0; i<10; ++i)
-    ks.add(100+10*i,true);
-
-  suite.err() << ks.score() << std::endl;
-
-  ks.p_value();
+  test_one_sample(suite);
+  test_p_value(suite);
 
   return suite.return_value();
 }
 
+void test_one_sample(test::Suite& suite)
+{
+  std::vector<double> correct(11);
+  for (size_t i=0; i<correct.size(); ++i) {
+    double s1 = 1.0 - i/10.0;
+    double s2 = 0.0-i/10.0;
+    if (std::abs(s1)>std::abs(s2))
+      correct[i] = s1;
+    else
+      correct[i] = s2;
+  }
+
+  for (size_t i=0; i<11; ++i) {
+    statistics::KolmogorovSmirnov ks;
+    for (size_t j=0; j<11; ++j) {
+      ks.add(j, i==j);
+    }     
+    double score = ks.signed_score();
+    if (!suite.add(suite.equal(score, correct[i]))) {
+      suite.err() << "signed_score(void) failed\n";
+    }
+  }
+
+  statistics::KolmogorovSmirnov ks;
+  for (size_t i=0; i<11; ++i) {
+    ks.add(i, i==0);
+  }     
+  size_t n=110000;
+  double p = ks.p_value(n);
+  double p_correct = 2.0/11.0;
+  double margin = 5*std::sqrt(n*p_correct);
+  if (p>p_correct+margin || p<p_correct-margin) {
+    suite.err() << "Error: p-value: " << p << "\n"
+                << "expected approximately: " << p_correct << "\n"
+                << "and at most " << margin << "deviation\n";
+    suite.add(false);
+  } 
+}
+
+void test_p_value(test::Suite& suite)
+{
+  statistics::KolmogorovSmirnov ks;
+  for (size_t i=0; i<100; ++i) {
+    ks.add(i, true);
+    ks.add(i+14.5, false);
+  } 
+  statistics::Averager a;
+  for (size_t n=0; n<100; ++n) {
+    a.add(ks.p_value(100));
+  }
+  double margin = 5 * a.standard_error(); 
+  double p_approx = ks.p_value();
+  if (std::abs(a.mean()-p_approx)>margin) {
+    suite.add(false);
+    suite.err() << "Error: unexpected large deviation between p_values\n"
+                << "permutation p-value: " << a.mean() << "\n"
+                << "analytical approximation: " << p_approx << "\n"
+                << "expected deviation to be smaller than " << margin << "\n";
+  }
+  
+}

trunk/yat/statistics/KolmogorovSmirnov.cc

@@ -61 +61 @@
     double res=0;
     double res2=0;
-    double s2 = score() * sqrt(sum_w1_*sum_w2_/(sum_w1_+sum_w2_));
-    s2 *= s2;
+    double s2 = score() * score() * sum_w1_*sum_w2_/(sum_w1_+sum_w2_);
     int sign = 1;
     for (size_t k = 1; k<100; ++k) {
-      res += sign * exp(-2.0 * k * k * s2);
+      res += sign * 2 * exp(-2.0 * k * k * s2);
       sign *= -1;
+      // ignore remaining terms as they are small
       if (res==res2)
         return res;

trunk/yat/statistics/KolmogorovSmirnov.h

@@ -99 +99 @@
     mutable double score_;
     typedef std::pair<std::pair<double,double>, bool> trip;
-    typedef std::set<trip, std::greater<trip> > data_w;
+    typedef std::multiset<trip> data_w;
     data_w data_;
     double sum_w1_;