source: trunk/yat/statistics/KolmogorovSmirnov.cc @ 1003

Last change on this file since 1003 was 1003, checked in by Peter, 14 years ago

adding class for Kolmogorov Smirnov test

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 2.4 KB
Line 
1// $Id: KolmogorovSmirnov.cc 1003 2008-01-18 22:13:13Z peter $
2
3#include "KolmogorovSmirnov.h"
4#include "yat/random/random.h"
5#include "yat/utility/stl_utility.h"
6
7#include <algorithm>
8#include <cassert>
9#include <cmath>
10#include <deque>
11#include <functional>
12
13
14namespace theplu {
15namespace yat {
16namespace statistics {
17
18
19  KolmogorovSmirnov::KolmogorovSmirnov(void)
20    : cached_(true), score_(0.0), sum_w1_(0.0), sum_w2_(0.0)
21  {
22  }
23
24
25  void KolmogorovSmirnov::add(double x, bool target, double w)
26  {
27    if (w==0) // ignoring zero weight data
28      return;
29
30    data_.insert(std::make_pair(std::make_pair(x,w),target));
31    if (target){
32      ++n1_;
33      sum_w1_+=w;
34    }
35    else {
36      ++n2_;
37      sum_w2_+=w;
38    }
39    cached_=false;
40  }
41
42
43  double KolmogorovSmirnov::p_value(size_t perm) const
44  {
45    size_t count=0;
46    std::deque<bool> targets;
47    typedef data_w::const_iterator const_iter;
48    for (const_iter i(data_.begin()); i!=data_.end(); ++i)
49      targets.push_back(i->second);
50
51    for (size_t i=0; i<perm; ++i){
52      yat::random::DiscreteUniform d;
53      std::random_shuffle(targets.begin(), targets.end(), d); 
54      KolmogorovSmirnov ks;
55      std::deque<bool>::const_iterator target_i(targets.begin());
56      const_iter end(data_.end());
57      for (const_iter i(data_.begin()); i!=end; ++i, ++target_i) 
58        ks.add(i->first.first, *target_i, i->first.second);
59     
60      if (ks.score()>=score())
61        ++count;
62    }
63    return static_cast<double>(count)/perm;
64  }
65
66
67  double KolmogorovSmirnov::score(void) const
68  {
69    if (cached_)
70      return score_;
71    if (!sum_w1_ || !sum_w2_)
72      return 0.0;
73    score_=0;
74    std::vector<double> v;
75    scores(v);
76    std::less<double> l;
77    utility::abs<double> a;
78    using utility::make_compose_f_gx_hy;
79    score_ = std::abs(*std::max_element(v.begin(), v.end(), 
80                                        make_compose_f_gx_hy(l, a, a)));
81    return score_;
82  }
83
84
85  void KolmogorovSmirnov::scores(std::vector<double>& res) const
86  {
87    res.clear();
88    res.reserve(data_.size());
89    data_w::const_iterator iter(data_.begin());
90    double f1=0;
91    double f2=0;
92    // Peter, we should take care of ties!
93    while(iter!=data_.end()){
94      if (iter->second)
95        f1 += iter->first.second;
96      else
97        f2 += iter->first.second;
98      res.push_back(f1/sum_w1_-f2/sum_w2_);
99    }
100  }
101
102
103  void KolmogorovSmirnov::reset(void)
104  {
105    cached_=false;
106    sum_w1_=0;
107    sum_w2_=0;
108  }
109
110
111  std::ostream& operator<<(std::ostream& os, const KolmogorovSmirnov& ks)
112  {
113    std::vector<double> s;
114    ks.scores(s);
115    for (size_t i=0; i<s.size(); ++i)
116      os << i << "\t" << s[i] << "\n";
117    return os;
118  }
119
120
121}}} // of namespace theplu yat statistics
Note: See TracBrowser for help on using the repository browser.