test/kolmogorov_smirnov.cc

Code
Comments
Other
Rev Date Author Line
1593 21 Oct 08 peter 1 // $Id$
1593 21 Oct 08 peter 2
1593 21 Oct 08 peter 3 /*
2998 14 Mar 13 peter 4   Copyright (C) 2008, 2009, 2010, 2012, 2013 Peter Johansson
1593 21 Oct 08 peter 5
1593 21 Oct 08 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
1593 21 Oct 08 peter 7
1593 21 Oct 08 peter 8   The yat library is free software; you can redistribute it and/or
1593 21 Oct 08 peter 9   modify it under the terms of the GNU General Public License as
1593 21 Oct 08 peter 10   published by the Free Software Foundation; either version 3 of the
1593 21 Oct 08 peter 11   License, or (at your option) any later version.
1593 21 Oct 08 peter 12
1593 21 Oct 08 peter 13   The yat library is distributed in the hope that it will be useful,
1593 21 Oct 08 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
1593 21 Oct 08 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1593 21 Oct 08 peter 16   General Public License for more details.
1593 21 Oct 08 peter 17
1593 21 Oct 08 peter 18   You should have received a copy of the GNU General Public License
1593 21 Oct 08 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
1593 21 Oct 08 peter 20 */
1593 21 Oct 08 peter 21
2881 18 Nov 12 peter 22 #include <config.h>
2881 18 Nov 12 peter 23
1593 21 Oct 08 peter 24 #include "Suite.h"
1593 21 Oct 08 peter 25
1701 07 Jan 09 peter 26 #include "yat/random/random.h"
1600 27 Oct 08 peter 27 #include "yat/statistics/Averager.h"
1593 21 Oct 08 peter 28 #include "yat/statistics/KolmogorovSmirnov.h"
2998 14 Mar 13 peter 29 #include "yat/statistics/KolmogorovSmirnovOneSample.h"
1593 21 Oct 08 peter 30
2202 21 Feb 10 peter 31 #include <boost/concept_archetype.hpp>
2202 21 Feb 10 peter 32
1600 27 Oct 08 peter 33 #include <cmath>
1701 07 Jan 09 peter 34 #include <deque>
1600 27 Oct 08 peter 35 #include <iostream>
1701 07 Jan 09 peter 36 #include <vector>
1593 21 Oct 08 peter 37
1600 27 Oct 08 peter 38 using namespace theplu::yat;
2721 12 Apr 12 peter 39 using statistics::KolmogorovSmirnov;
1600 27 Oct 08 peter 40
1600 27 Oct 08 peter 41 void test_one_sample(test::Suite&);
1617 06 Nov 08 peter 42 void test_two_sample(test::Suite&);
1600 27 Oct 08 peter 43 void test_p_value(test::Suite&);
1701 07 Jan 09 peter 44 void test_shuffle(test::Suite&);
1701 07 Jan 09 peter 45 void test_range(test::Suite&);
2721 12 Apr 12 peter 46 void test_remove(test::Suite&);
1612 04 Nov 08 peter 47 void test_reset(test::Suite&);
1687 30 Dec 08 peter 48 void test_ties(test::Suite&);
2202 21 Feb 10 peter 49 void test_compile(void);
2998 14 Mar 13 peter 50 void test_ks_one_sample(test::Suite&);
1600 27 Oct 08 peter 51
1593 21 Oct 08 peter 52 int main(int argc, char* argv[])
2721 12 Apr 12 peter 53 {
1593 21 Oct 08 peter 54   test::Suite suite(argc, argv);
1593 21 Oct 08 peter 55
1600 27 Oct 08 peter 56   test_one_sample(suite);
1617 06 Nov 08 peter 57   test_two_sample(suite);
1626 15 Nov 08 peter 58   test_p_value(suite);
1701 07 Jan 09 peter 59   test_shuffle(suite);
1701 07 Jan 09 peter 60   test_range(suite);
1612 04 Nov 08 peter 61   test_reset(suite);
1687 30 Dec 08 peter 62   test_ties(suite);
2202 21 Feb 10 peter 63   test_compile();
2721 12 Apr 12 peter 64   test_remove(suite);
2998 14 Mar 13 peter 65   test_ks_one_sample(suite);
1593 21 Oct 08 peter 66
1600 27 Oct 08 peter 67   return suite.return_value();
1600 27 Oct 08 peter 68 }
1593 21 Oct 08 peter 69
2202 21 Feb 10 peter 70 void test_compile(void)
2202 21 Feb 10 peter 71 {
2202 21 Feb 10 peter 72   using statistics::KolmogorovSmirnov;
2202 21 Feb 10 peter 73   // do not run compiler test
2202 21 Feb 10 peter 74   if (false) {
2202 21 Feb 10 peter 75     KolmogorovSmirnov ks;
2202 21 Feb 10 peter 76     ks.add(boost::forward_iterator_archetype<KolmogorovSmirnov::Element>(),
2202 21 Feb 10 peter 77            boost::forward_iterator_archetype<KolmogorovSmirnov::Element>());
2202 21 Feb 10 peter 78   }
2202 21 Feb 10 peter 79 }
2202 21 Feb 10 peter 80
1600 27 Oct 08 peter 81 void test_one_sample(test::Suite& suite)
1600 27 Oct 08 peter 82 {
1600 27 Oct 08 peter 83   std::vector<double> correct(11);
1600 27 Oct 08 peter 84   for (size_t i=0; i<correct.size(); ++i) {
1600 27 Oct 08 peter 85     double s1 = 1.0 - i/10.0;
1600 27 Oct 08 peter 86     double s2 = 0.0-i/10.0;
1600 27 Oct 08 peter 87     if (std::abs(s1)>std::abs(s2))
1600 27 Oct 08 peter 88       correct[i] = s1;
1600 27 Oct 08 peter 89     else
1600 27 Oct 08 peter 90       correct[i] = s2;
1600 27 Oct 08 peter 91   }
1593 21 Oct 08 peter 92
1600 27 Oct 08 peter 93   for (size_t i=0; i<11; ++i) {
1600 27 Oct 08 peter 94     statistics::KolmogorovSmirnov ks;
1600 27 Oct 08 peter 95     for (size_t j=0; j<11; ++j) {
1600 27 Oct 08 peter 96       ks.add(j, i==j);
2721 12 Apr 12 peter 97     }
1600 27 Oct 08 peter 98     double score = ks.signed_score();
1600 27 Oct 08 peter 99     if (!suite.add(suite.equal(score, correct[i]))) {
1600 27 Oct 08 peter 100       suite.err() << "signed_score(void) failed\n";
1600 27 Oct 08 peter 101     }
1600 27 Oct 08 peter 102   }
1600 27 Oct 08 peter 103
1600 27 Oct 08 peter 104   statistics::KolmogorovSmirnov ks;
1600 27 Oct 08 peter 105   for (size_t i=0; i<11; ++i) {
1600 27 Oct 08 peter 106     ks.add(i, i==0);
2721 12 Apr 12 peter 107   }
1600 27 Oct 08 peter 108   size_t n=110000;
1600 27 Oct 08 peter 109   double p = ks.p_value(n);
1600 27 Oct 08 peter 110   double p_correct = 2.0/11.0;
1618 06 Nov 08 peter 111   double margin = 10*std::sqrt(p_correct*(1-p_correct)/n);
1600 27 Oct 08 peter 112   if (p>p_correct+margin || p<p_correct-margin) {
1600 27 Oct 08 peter 113     suite.err() << "Error: p-value: " << p << "\n"
1600 27 Oct 08 peter 114                 << "expected approximately: " << p_correct << "\n"
1600 27 Oct 08 peter 115                 << "and at most " << margin << "deviation\n";
1600 27 Oct 08 peter 116     suite.add(false);
2721 12 Apr 12 peter 117   }
1593 21 Oct 08 peter 118 }
1593 21 Oct 08 peter 119
1617 06 Nov 08 peter 120
1617 06 Nov 08 peter 121 void test_two_sample(test::Suite& suite)
1617 06 Nov 08 peter 122 {
1617 06 Nov 08 peter 123   suite.err() << "testing two sample\n";
1617 06 Nov 08 peter 124   statistics::KolmogorovSmirnov ks;
1617 06 Nov 08 peter 125   for (size_t i=0; i<5; ++i)
1617 06 Nov 08 peter 126     ks.add(i+0.5, i<2);
1617 06 Nov 08 peter 127   suite.add(suite.equal(ks.score(), 1.0));
1618 06 Nov 08 peter 128   size_t n=100000;
1617 06 Nov 08 peter 129   double p = ks.p_value(n);
1617 06 Nov 08 peter 130   double p_correct = 0.2;
1618 06 Nov 08 peter 131   double margin=10*std::sqrt(p_correct*(1-p_correct)/n);
1704 08 Jan 09 peter 132   if (!suite.equal_fix(p, p_correct, margin) ) {
1617 06 Nov 08 peter 133     suite.add(false);
1617 06 Nov 08 peter 134     suite.err() << "Error: p = " << p << "\n"
1617 06 Nov 08 peter 135                 << "correct p would be: " << p_correct << "\n"
1617 06 Nov 08 peter 136                 << "expected a difference less than " << margin << "margin\n";
1617 06 Nov 08 peter 137     suite.err() << p << std::endl;
1617 06 Nov 08 peter 138   }
1617 06 Nov 08 peter 139 }
1617 06 Nov 08 peter 140
1617 06 Nov 08 peter 141
1600 27 Oct 08 peter 142 void test_p_value(test::Suite& suite)
1600 27 Oct 08 peter 143 {
1600 27 Oct 08 peter 144   statistics::KolmogorovSmirnov ks;
1600 27 Oct 08 peter 145   for (size_t i=0; i<100; ++i) {
1600 27 Oct 08 peter 146     ks.add(i, true);
1600 27 Oct 08 peter 147     ks.add(i+14.5, false);
2721 12 Apr 12 peter 148   }
1626 15 Nov 08 peter 149   suite.add(suite.equal(ks.score(), 0.15, 10));
1626 15 Nov 08 peter 150
1600 27 Oct 08 peter 151   statistics::Averager a;
1600 27 Oct 08 peter 152   for (size_t n=0; n<100; ++n) {
1600 27 Oct 08 peter 153     a.add(ks.p_value(100));
1600 27 Oct 08 peter 154   }
2721 12 Apr 12 peter 155   double margin = 5 * a.standard_error();
1600 27 Oct 08 peter 156   double p_approx = ks.p_value();
1704 08 Jan 09 peter 157   if (!suite.equal_fix(a.mean(), p_approx, margin) ) {
1600 27 Oct 08 peter 158     suite.add(false);
1600 27 Oct 08 peter 159     suite.err() << "Error: unexpected large deviation between p_values\n"
1600 27 Oct 08 peter 160                 << "permutation p-value: " << a.mean() << "\n"
1600 27 Oct 08 peter 161                 << "analytical approximation: " << p_approx << "\n"
1600 27 Oct 08 peter 162                 << "expected deviation to be smaller than " << margin << "\n";
1600 27 Oct 08 peter 163   }
2721 12 Apr 12 peter 164
1600 27 Oct 08 peter 165 }
1612 04 Nov 08 peter 166
1701 07 Jan 09 peter 167
1701 07 Jan 09 peter 168 void test_range(test::Suite& suite)
1701 07 Jan 09 peter 169 {
1701 07 Jan 09 peter 170   suite.err() << "testing range" << std::endl;
1701 07 Jan 09 peter 171   std::deque<bool> labels;
1701 07 Jan 09 peter 172   for (size_t i=0; i<10; ++i) {
1701 07 Jan 09 peter 173     labels.push_back(i<5);
1701 07 Jan 09 peter 174   }
1701 07 Jan 09 peter 175   std::vector<statistics::KolmogorovSmirnov::Element> data;
1701 07 Jan 09 peter 176   statistics::KolmogorovSmirnov::Element elem;
1701 07 Jan 09 peter 177   elem.weight = 1.0;
1701 07 Jan 09 peter 178   for (size_t i=0; i<10; ++i) {
1701 07 Jan 09 peter 179     elem.value = i;
1701 07 Jan 09 peter 180     elem.label = labels[i];
1701 07 Jan 09 peter 181     data.push_back(elem);
1701 07 Jan 09 peter 182   }
1701 07 Jan 09 peter 183   statistics::KolmogorovSmirnov ks;
1701 07 Jan 09 peter 184   ks.add(data.begin(), data.end());
1701 07 Jan 09 peter 185   suite.add(suite.equal(ks.score(), 1.0));
2721 12 Apr 12 peter 186
1701 07 Jan 09 peter 187   // testing that adding a range gives same result as adding elements
1701 07 Jan 09 peter 188   // sequentially
1701 07 Jan 09 peter 189   statistics::KolmogorovSmirnov ks2;
1701 07 Jan 09 peter 190   for (size_t i=0; i<data.size(); ++i)
1701 07 Jan 09 peter 191     ks2.add(data[i].value, data[i].label, data[i].weight);
1701 07 Jan 09 peter 192   suite.add(suite.equal(ks2.signed_score(), ks.signed_score()));
1701 07 Jan 09 peter 193
1701 07 Jan 09 peter 194   theplu::yat::random::random_shuffle(labels.begin(), labels.end());
1701 07 Jan 09 peter 195   for (size_t i=0; i<data.size(); ++i) {
1701 07 Jan 09 peter 196     data[i].label = labels[i];
1701 07 Jan 09 peter 197   }
1701 07 Jan 09 peter 198   ks.reset();
1701 07 Jan 09 peter 199   ks.add(data.begin(), data.end());
1701 07 Jan 09 peter 200   ks2.reset();
1701 07 Jan 09 peter 201   for (size_t i=0; i<data.size(); ++i)
1701 07 Jan 09 peter 202     ks2.add(data[i].value, data[i].label, data[i].weight);
1701 07 Jan 09 peter 203   suite.add(suite.equal(ks2.signed_score(), ks.signed_score()));
1701 07 Jan 09 peter 204 }
1701 07 Jan 09 peter 205
1701 07 Jan 09 peter 206
2721 12 Apr 12 peter 207 void test_remove(test::Suite& suite)
2721 12 Apr 12 peter 208 {
2721 12 Apr 12 peter 209   suite.out() << "test remove\n";
2721 12 Apr 12 peter 210   KolmogorovSmirnov ks;
3018 04 Apr 13 peter 211   ks.add(0, true);
2721 12 Apr 12 peter 212   ks.add(1, true);
2721 12 Apr 12 peter 213   ks.add(2, false);
2721 12 Apr 12 peter 214   ks.add(2, true);
2721 12 Apr 12 peter 215   ks.add(3, false);
2721 12 Apr 12 peter 216   double score = ks.score();
3018 04 Apr 13 peter 217   double x = 0;
3018 04 Apr 13 peter 218   ks.add(x, false);
2721 12 Apr 12 peter 219
2721 12 Apr 12 peter 220   try {
3018 04 Apr 13 peter 221     ks.remove(x, false);
2721 12 Apr 12 peter 222   }
2721 12 Apr 12 peter 223   catch (std::runtime_error& e) {
2721 12 Apr 12 peter 224     suite.out() << "what(): " << e.what() << "\n";
2721 12 Apr 12 peter 225     suite.add(false);
2721 12 Apr 12 peter 226   }
2721 12 Apr 12 peter 227   suite.add(suite.equal(score, ks.score()));
2721 12 Apr 12 peter 228
2721 12 Apr 12 peter 229   try {
2721 12 Apr 12 peter 230     ks.remove(1,false);
2721 12 Apr 12 peter 231     suite.add(false);
2721 12 Apr 12 peter 232     suite.out() << "error: missing exception\n";
2721 12 Apr 12 peter 233   }
2721 12 Apr 12 peter 234   catch (std::runtime_error& e) {
2721 12 Apr 12 peter 235     suite.add(true);
2721 12 Apr 12 peter 236   }
2721 12 Apr 12 peter 237 }
2721 12 Apr 12 peter 238
2721 12 Apr 12 peter 239
1612 04 Nov 08 peter 240 void test_reset(test::Suite& suite)
1612 04 Nov 08 peter 241 {
1612 04 Nov 08 peter 242   suite.err() << "testing reset\n";
1612 04 Nov 08 peter 243   statistics::KolmogorovSmirnov ks;
1612 04 Nov 08 peter 244   ks.add(1.0, true);
1612 04 Nov 08 peter 245   ks.add(2.0, false);
1612 04 Nov 08 peter 246   ks.add(3.0, true);
1612 04 Nov 08 peter 247   double score = ks.score();
1612 04 Nov 08 peter 248   double p = ks.p_value();
1612 04 Nov 08 peter 249   ks.reset();
1612 04 Nov 08 peter 250   ks.add(1.0, true);
1612 04 Nov 08 peter 251   ks.add(2.0, false);
1612 04 Nov 08 peter 252   ks.add(3.0, true);
1612 04 Nov 08 peter 253   suite.add(suite.equal(ks.score(), score));
1612 04 Nov 08 peter 254   suite.add(suite.equal(ks.p_value(), p));
1612 04 Nov 08 peter 255 }
1687 30 Dec 08 peter 256
1687 30 Dec 08 peter 257
1687 30 Dec 08 peter 258 void test_ties(test::Suite& suite)
1687 30 Dec 08 peter 259 {
1687 30 Dec 08 peter 260   suite.err() << "test ties" << std::endl;
1687 30 Dec 08 peter 261   statistics::KolmogorovSmirnov ks;
1687 30 Dec 08 peter 262   for (size_t i=0; i<5; ++i)
1687 30 Dec 08 peter 263     ks.add(i, true);
1687 30 Dec 08 peter 264   ks.add(0, false);
1688 30 Dec 08 peter 265   suite.add(suite.equal(ks.score(), 1.0-0.2));
1687 30 Dec 08 peter 266 }
1701 07 Jan 09 peter 267
1701 07 Jan 09 peter 268
1701 07 Jan 09 peter 269 void test_shuffle(test::Suite& suite)
1701 07 Jan 09 peter 270 {
1701 07 Jan 09 peter 271   suite.err() << "testing shuffle" << std::endl;
1701 07 Jan 09 peter 272   statistics::KolmogorovSmirnov ks;
1701 07 Jan 09 peter 273   for (size_t i=0; i<10; ++i)
1701 07 Jan 09 peter 274     ks.add(i, i<5);
1701 07 Jan 09 peter 275   ks.shuffle();
1701 07 Jan 09 peter 276 }
2998 14 Mar 13 peter 277
2998 14 Mar 13 peter 278
2998 14 Mar 13 peter 279 void test_ks_one_sample(test::Suite& suite)
2998 14 Mar 13 peter 280 {
2998 14 Mar 13 peter 281   suite.err() << "testing one sample" << std::endl;
2998 14 Mar 13 peter 282   statistics::KolmogorovSmirnovOneSample ks;
2998 14 Mar 13 peter 283   for (size_t i=0; i<=10; ++i)
2998 14 Mar 13 peter 284     ks.add(i/10.0);
2998 14 Mar 13 peter 285   suite.out() << "score: " << ks.score() << "\n";
2998 14 Mar 13 peter 286   suite.out() << "P: " << ks.p_value() << "\n";
2998 14 Mar 13 peter 287   ks.remove(0);
2998 14 Mar 13 peter 288   ks.remove(0.1);
2998 14 Mar 13 peter 289   ks.remove(0.2);
2998 14 Mar 13 peter 290   ks.remove(0.3);
2998 14 Mar 13 peter 291   suite.out() << "P: " << ks.p_value() << "\n";
2998 14 Mar 13 peter 292   suite.out() << "score: " << ks.score() << "\n";
2998 14 Mar 13 peter 293   suite.out() << "score: " << ks.signed_score() << "\n";
2998 14 Mar 13 peter 294   ks.reset();
2998 14 Mar 13 peter 295 }

Generated on Sat Oct 14 02:21:06 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)