source: trunk/test/kolmogorov_smirnov.cc @ 3018

// $Id: kolmogorov_smirnov.cc 3018 2013-04-04 04:46:38Z peter $

/*
  Copyright (C) 2008, 2009, 2010, 2012, 2013 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 <config.h>

#include "Suite.h"

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

#include <boost/concept_archetype.hpp>

#include <cmath>
#include <deque>
#include <iostream>
#include <vector>

using namespace theplu::yat;
using statistics::KolmogorovSmirnov;

void test_one_sample(test::Suite&);
void test_two_sample(test::Suite&);
void test_p_value(test::Suite&);
void test_shuffle(test::Suite&);
void test_range(test::Suite&);
void test_remove(test::Suite&);
void test_reset(test::Suite&);
void test_ties(test::Suite&);
void test_compile(void);
void test_ks_one_sample(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_shuffle(suite);
  test_range(suite);
  test_reset(suite);
  test_ties(suite);
  test_compile();
  test_remove(suite);
  test_ks_one_sample(suite);

  return suite.return_value();
}

void test_compile(void)
{
  using statistics::KolmogorovSmirnov;
  // do not run compiler test
  if (false) {
    KolmogorovSmirnov ks;
    ks.add(boost::forward_iterator_archetype<KolmogorovSmirnov::Element>(),
           boost::forward_iterator_archetype<KolmogorovSmirnov::Element>());
  }
}

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 (!suite.equal_fix(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 (!suite.equal_fix(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_range(test::Suite& suite)
{
  suite.err() << "testing range" << std::endl;
  std::deque<bool> labels;
  for (size_t i=0; i<10; ++i) {
    labels.push_back(i<5);
  }
  std::vector<statistics::KolmogorovSmirnov::Element> data;
  statistics::KolmogorovSmirnov::Element elem;
  elem.weight = 1.0;
  for (size_t i=0; i<10; ++i) {
    elem.value = i;
    elem.label = labels[i];
    data.push_back(elem);
  }
  statistics::KolmogorovSmirnov ks;
  ks.add(data.begin(), data.end());
  suite.add(suite.equal(ks.score(), 1.0));

  // testing that adding a range gives same result as adding elements
  // sequentially
  statistics::KolmogorovSmirnov ks2;
  for (size_t i=0; i<data.size(); ++i)
    ks2.add(data[i].value, data[i].label, data[i].weight);
  suite.add(suite.equal(ks2.signed_score(), ks.signed_score()));

  theplu::yat::random::random_shuffle(labels.begin(), labels.end());
  for (size_t i=0; i<data.size(); ++i) {
    data[i].label = labels[i];
  }
  ks.reset();
  ks.add(data.begin(), data.end());
  ks2.reset();
  for (size_t i=0; i<data.size(); ++i)
    ks2.add(data[i].value, data[i].label, data[i].weight);
  suite.add(suite.equal(ks2.signed_score(), ks.signed_score()));
}


void test_remove(test::Suite& suite)
{
  suite.out() << "test remove\n";
  KolmogorovSmirnov ks;
  ks.add(0, true);
  ks.add(1, true);
  ks.add(2, false);
  ks.add(2, true);
  ks.add(3, false);
  double score = ks.score();
  double x = 0;
  ks.add(x, false);

  try {
    ks.remove(x, false);
  }
  catch (std::runtime_error& e) {
    suite.out() << "what(): " << e.what() << "\n";
    suite.add(false);
  }
  suite.add(suite.equal(score, ks.score()));

  try {
    ks.remove(1,false);
    suite.add(false);
    suite.out() << "error: missing exception\n";
  }
  catch (std::runtime_error& e) {
    suite.add(true);
  }
}


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


void test_ties(test::Suite& suite)
{
  suite.err() << "test ties" << std::endl;
  statistics::KolmogorovSmirnov ks;
  for (size_t i=0; i<5; ++i)
    ks.add(i, true);
  ks.add(0, false);
  suite.add(suite.equal(ks.score(), 1.0-0.2));
}


void test_shuffle(test::Suite& suite)
{
  suite.err() << "testing shuffle" << std::endl;
  statistics::KolmogorovSmirnov ks;
  for (size_t i=0; i<10; ++i)
    ks.add(i, i<5);
  ks.shuffle();
}


void test_ks_one_sample(test::Suite& suite)
{
  suite.err() << "testing one sample" << std::endl;
  statistics::KolmogorovSmirnovOneSample ks;
  for (size_t i=0; i<=10; ++i)
    ks.add(i/10.0);
  suite.out() << "score: " << ks.score() << "\n";
  suite.out() << "P: " << ks.p_value() << "\n";
  ks.remove(0);
  ks.remove(0.1);
  ks.remove(0.2);
  ks.remove(0.3);
  suite.out() << "P: " << ks.p_value() << "\n";
  suite.out() << "score: " << ks.score() << "\n";
  suite.out() << "score: " << ks.signed_score() << "\n";
  ks.reset();
}