source: trunk/test/kolmogorov_smirnov_test.cc @ 1626

// $Id: kolmogorov_smirnov_test.cc 1626 2008-11-15 16:01:46Z peter $

/*
  Copyright (C) 2008 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "Suite.h"

#include "yat/statistics/Averager.h"
#include "yat/statistics/KolmogorovSmirnov.h"
#include "yat/random/random.h"

#include <cmath>
#include <iostream>

using namespace theplu::yat;

void test_one_sample(test::Suite&);
void test_two_sample(test::Suite&);
void test_p_value(test::Suite&);
void test_reset(test::Suite&);

int main(int argc, char* argv[])
{  
  test::Suite suite(argc, argv);

  test_one_sample(suite);
  test_two_sample(suite);
  test_p_value(suite);
  test_reset(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 = 10*std::sqrt(p_correct*(1-p_correct)/n);
  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_two_sample(test::Suite& suite)
{
  suite.err() << "testing two sample\n";
  statistics::KolmogorovSmirnov ks;
  for (size_t i=0; i<5; ++i)
    ks.add(i+0.5, i<2);
  suite.add(suite.equal(ks.score(), 1.0));
  size_t n=100000;
  double p = ks.p_value(n);
  double p_correct = 0.2;
  double margin=10*std::sqrt(p_correct*(1-p_correct)/n);
  if (std::abs(p-p_correct)>margin) {
    suite.add(false);
    suite.err() << "Error: p = " << p << "\n"
                << "correct p would be: " << p_correct << "\n"
                << "expected a difference less than " << margin << "margin\n";
    suite.err() << p << std::endl;
  }
}


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);
  } 
  suite.add(suite.equal(ks.score(), 0.15, 10));

  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";
  }
  
}

void test_reset(test::Suite& suite)
{
  suite.err() << "testing reset\n";
  statistics::KolmogorovSmirnov ks;
  ks.add(1.0, true);
  ks.add(2.0, false);
  ks.add(3.0, true);
  double score = ks.score();
  double p = ks.p_value();
  ks.reset();
  ks.add(1.0, true);
  ks.add(2.0, false);
  ks.add(3.0, true);
  suite.add(suite.equal(ks.score(), score));
  suite.add(suite.equal(ks.p_value(), p));
}