source: trunk/test/kolmogorov_smirnov.cc @ 2881

Last change on this file since 2881 was 2881, checked in by Peter, 9 years ago

Define PP variables in config.h rather than in CPPFLAGS. Include
config.h into all source files. Only ammend CXXFLAGS with '-Wall
-pedantic' when --enable-debug. In default mode we respect CXXFLAGS
value set by user, or set to default value '-O3'.

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