source: trunk/yat/statistics/KolmogorovSmirnovOneSample.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: 2.9 KB
Line 
1// $Id: KolmogorovSmirnovOneSample.cc 2998 2013-03-14 08:08:54Z peter $
2
3/*
4  Copyright (C) 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 "KolmogorovSmirnovOneSample.h"
25
26#include "yat/utility/Exception.h"
27
28#include <cassert>
29#include <cmath>
30#include <sstream>
31
32#include <iostream>
33
34namespace theplu {
35namespace yat {
36namespace statistics {
37
38  KolmogorovSmirnovOneSample::KolmogorovSmirnovOneSample(void)
39    : sum_w_(0)
40  {
41  }
42
43
44  void KolmogorovSmirnovOneSample::add(double value, double weight)
45  {
46    if (weight==0)
47      return;
48    cached_ = false;
49    data_.insert(std::make_pair(value, weight));
50    sum_w_ += weight;
51  }
52
53
54  double KolmogorovSmirnovOneSample::p_value(void) const
55  {
56    double res=0;
57    double res2=0;
58    double s2 = score() * score() * sum_w_;
59    int sign = 1;
60    for (size_t k = 1; k<100; ++k) {
61      res += sign * 2 * exp(-2.0 * k * k * s2);
62      sign *= -1;
63      // ignore remaining terms as they are small
64      if (res==res2)
65        return res;
66      res2 = res;
67    }
68
69    return res;
70  }
71
72
73  void KolmogorovSmirnovOneSample::remove(double value, double weight)
74  {
75    if (!weight)
76      return;
77    typedef std::multimap<double, double>::iterator iterator;
78    std::pair<iterator, iterator> range = data_.equal_range(value);
79    for ( ; range.first!=range.second; ++range.first) {
80      if (range.first->second == weight) {
81        data_.erase(range.first);
82        sum_w_ -= weight;
83        cached_ = false;
84        return;
85      }
86    }
87    std::ostringstream ss;
88    ss << "KolmogorovSmirnovOneSample::remove: " << value << " " << weight;
89    throw utility::runtime_error(ss.str());
90  }
91
92
93  void KolmogorovSmirnovOneSample::reset(void)
94  {
95    data_.clear();
96    sum_w_ = 0;
97  }
98
99
100  double KolmogorovSmirnovOneSample::score(void) const
101  {
102    return std::abs(signed_score());
103  }
104
105
106  double KolmogorovSmirnovOneSample::signed_score(void) const
107  {
108    if (cached_)
109      return score_;
110    score_ = 0.0;
111    if (data_.empty())
112      return 0.0;
113
114    typedef std::multimap<double, double>::const_iterator iterator;
115    iterator iter = data_.begin();
116    double w = 0;
117    while (iter != data_.end()) {
118      double x = iter->first;
119      while (iter!=data_.end() && iter->first==x) {
120        w += iter->second;
121        ++iter;
122      }
123      double s = w / sum_w_ - x;
124      if (std::abs(s) > std::abs(score_))
125        score_ = s;
126    }
127
128    cached_ = true;
129    return score_;
130  }
131
132}}} // of namespace theplu yat statistics
Note: See TracBrowser for help on using the repository browser.