source: trunk/test/kolmogorov_smirnov.cc @ 2998

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

New class KolmogorovSmirnovOneSample?. closes #686

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.2 KB
Line 
1// $Id: kolmogorov_smirnov.cc 2998 2013-03-14 08:08:54Z peter $
2
3/*
4  Copyright (C) 2008, 2009, 2010, 2012, 2013 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#include "yat/statistics/KolmogorovSmirnovOneSample.h"
30
31#include <boost/concept_archetype.hpp>
32
33#include <cmath>
34#include <deque>
35#include <iostream>
36#include <vector>
37
38using namespace theplu::yat;
39using statistics::KolmogorovSmirnov;
40
41void test_one_sample(test::Suite&);
42void test_two_sample(test::Suite&);
43void test_p_value(test::Suite&);
44void test_shuffle(test::Suite&);
45void test_range(test::Suite&);
46void test_remove(test::Suite&);
47void test_reset(test::Suite&);
48void test_ties(test::Suite&);
49void test_compile(void);
50void test_ks_one_sample(test::Suite&);
51
52int main(int argc, char* argv[])
53{
54  test::Suite suite(argc, argv);
55
56  test_one_sample(suite);
57  test_two_sample(suite);
58  test_p_value(suite);
59  test_shuffle(suite);
60  test_range(suite);
61  test_reset(suite);
62  test_ties(suite);
63  test_compile();
64  test_remove(suite);
65  test_ks_one_sample(suite);
66
67  return suite.return_value();
68}
69
70void test_compile(void)
71{
72  using statistics::KolmogorovSmirnov;
73  // do not run compiler test
74  if (false) {
75    KolmogorovSmirnov ks;
76    ks.add(boost::forward_iterator_archetype<KolmogorovSmirnov::Element>(),
77           boost::forward_iterator_archetype<KolmogorovSmirnov::Element>());
78  }
79}
80
81void test_one_sample(test::Suite& suite)
82{
83  std::vector<double> correct(11);
84  for (size_t i=0; i<correct.size(); ++i) {
85    double s1 = 1.0 - i/10.0;
86    double s2 = 0.0-i/10.0;
87    if (std::abs(s1)>std::abs(s2))
88      correct[i] = s1;
89    else
90      correct[i] = s2;
91  }
92
93  for (size_t i=0; i<11; ++i) {
94    statistics::KolmogorovSmirnov ks;
95    for (size_t j=0; j<11; ++j) {
96      ks.add(j, i==j);
97    }
98    double score = ks.signed_score();
99    if (!suite.add(suite.equal(score, correct[i]))) {
100      suite.err() << "signed_score(void) failed\n";
101    }
102  }
103
104  statistics::KolmogorovSmirnov ks;
105  for (size_t i=0; i<11; ++i) {
106    ks.add(i, i==0);
107  }
108  size_t n=110000;
109  double p = ks.p_value(n);
110  double p_correct = 2.0/11.0;
111  double margin = 10*std::sqrt(p_correct*(1-p_correct)/n);
112  if (p>p_correct+margin || p<p_correct-margin) {
113    suite.err() << "Error: p-value: " << p << "\n"
114                << "expected approximately: " << p_correct << "\n"
115                << "and at most " << margin << "deviation\n";
116    suite.add(false);
117  }
118}
119
120
121void test_two_sample(test::Suite& suite)
122{
123  suite.err() << "testing two sample\n";
124  statistics::KolmogorovSmirnov ks;
125  for (size_t i=0; i<5; ++i)
126    ks.add(i+0.5, i<2);
127  suite.add(suite.equal(ks.score(), 1.0));
128  size_t n=100000;
129  double p = ks.p_value(n);
130  double p_correct = 0.2;
131  double margin=10*std::sqrt(p_correct*(1-p_correct)/n);
132  if (!suite.equal_fix(p, p_correct, margin) ) {
133    suite.add(false);
134    suite.err() << "Error: p = " << p << "\n"
135                << "correct p would be: " << p_correct << "\n"
136                << "expected a difference less than " << margin << "margin\n";
137    suite.err() << p << std::endl;
138  }
139}
140
141
142void test_p_value(test::Suite& suite)
143{
144  statistics::KolmogorovSmirnov ks;
145  for (size_t i=0; i<100; ++i) {
146    ks.add(i, true);
147    ks.add(i+14.5, false);
148  }
149  suite.add(suite.equal(ks.score(), 0.15, 10));
150
151  statistics::Averager a;
152  for (size_t n=0; n<100; ++n) {
153    a.add(ks.p_value(100));
154  }
155  double margin = 5 * a.standard_error();
156  double p_approx = ks.p_value();
157  if (!suite.equal_fix(a.mean(), p_approx, margin) ) {
158    suite.add(false);
159    suite.err() << "Error: unexpected large deviation between p_values\n"
160                << "permutation p-value: " << a.mean() << "\n"
161                << "analytical approximation: " << p_approx << "\n"
162                << "expected deviation to be smaller than " << margin << "\n";
163  }
164
165}
166
167
168void test_range(test::Suite& suite)
169{
170  suite.err() << "testing range" << std::endl;
171  std::deque<bool> labels;
172  for (size_t i=0; i<10; ++i) {
173    labels.push_back(i<5);
174  }
175  std::vector<statistics::KolmogorovSmirnov::Element> data;
176  statistics::KolmogorovSmirnov::Element elem;
177  elem.weight = 1.0;
178  for (size_t i=0; i<10; ++i) {
179    elem.value = i;
180    elem.label = labels[i];
181    data.push_back(elem);
182  }
183  statistics::KolmogorovSmirnov ks;
184  ks.add(data.begin(), data.end());
185  suite.add(suite.equal(ks.score(), 1.0));
186
187  // testing that adding a range gives same result as adding elements
188  // sequentially
189  statistics::KolmogorovSmirnov ks2;
190  for (size_t i=0; i<data.size(); ++i)
191    ks2.add(data[i].value, data[i].label, data[i].weight);
192  suite.add(suite.equal(ks2.signed_score(), ks.signed_score()));
193
194  theplu::yat::random::random_shuffle(labels.begin(), labels.end());
195  for (size_t i=0; i<data.size(); ++i) {
196    data[i].label = labels[i];
197  }
198  ks.reset();
199  ks.add(data.begin(), data.end());
200  ks2.reset();
201  for (size_t i=0; i<data.size(); ++i)
202    ks2.add(data[i].value, data[i].label, data[i].weight);
203  suite.add(suite.equal(ks2.signed_score(), ks.signed_score()));
204}
205
206
207void test_remove(test::Suite& suite)
208{
209  suite.out() << "test remove\n";
210  KolmogorovSmirnov ks;
211  ks.add(1, true);
212  ks.add(2, false);
213  ks.add(2, true);
214  ks.add(3, false);
215  double score = ks.score();
216  ks.add(4, false);
217
218  try {
219    ks.remove(4, false);
220  }
221  catch (std::runtime_error& e) {
222    suite.out() << "what(): " << e.what() << "\n";
223    suite.add(false);
224  }
225  suite.add(suite.equal(score, ks.score()));
226
227  try {
228    ks.remove(1,false);
229    suite.add(false);
230    suite.out() << "error: missing exception\n";
231  }
232  catch (std::runtime_error& e) {
233    suite.add(true);
234  }
235}
236
237
238void test_reset(test::Suite& suite)
239{
240  suite.err() << "testing reset\n";
241  statistics::KolmogorovSmirnov ks;
242  ks.add(1.0, true);
243  ks.add(2.0, false);
244  ks.add(3.0, true);
245  double score = ks.score();
246  double p = ks.p_value();
247  ks.reset();
248  ks.add(1.0, true);
249  ks.add(2.0, false);
250  ks.add(3.0, true);
251  suite.add(suite.equal(ks.score(), score));
252  suite.add(suite.equal(ks.p_value(), p));
253}
254
255
256void test_ties(test::Suite& suite)
257{
258  suite.err() << "test ties" << std::endl;
259  statistics::KolmogorovSmirnov ks;
260  for (size_t i=0; i<5; ++i)
261    ks.add(i, true);
262  ks.add(0, false);
263  suite.add(suite.equal(ks.score(), 1.0-0.2));
264}
265
266
267void test_shuffle(test::Suite& suite)
268{
269  suite.err() << "testing shuffle" << std::endl;
270  statistics::KolmogorovSmirnov ks;
271  for (size_t i=0; i<10; ++i)
272    ks.add(i, i<5);
273  ks.shuffle();
274}
275
276
277void test_ks_one_sample(test::Suite& suite)
278{
279  suite.err() << "testing one sample" << std::endl;
280  statistics::KolmogorovSmirnovOneSample ks;
281  for (size_t i=0; i<=10; ++i)
282    ks.add(i/10.0);
283  suite.out() << "score: " << ks.score() << "\n";
284  suite.out() << "P: " << ks.p_value() << "\n";
285  ks.remove(0);
286  ks.remove(0.1);
287  ks.remove(0.2);
288  ks.remove(0.3);
289  suite.out() << "P: " << ks.p_value() << "\n";
290  suite.out() << "score: " << ks.score() << "\n";
291  suite.out() << "score: " << ks.signed_score() << "\n";
292  ks.reset();
293}
Note: See TracBrowser for help on using the repository browser.