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

// $Id: KolmogorovSmirnov.cc 3018 2013-04-04 04:46:38Z peter $

/*
  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2011, 2012, 2013 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "KolmogorovSmirnov.h"
#include "yat/random/random.h"
#include "yat/utility/Exception.h"
#include "yat/utility/stl_utility.h"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <deque>
#include <functional>
#include <limits>
#include <sstream>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {


  KolmogorovSmirnov::Element::Element(void)
  {}


  KolmogorovSmirnov::Element::Element(double x, bool class_label, double w)
    : value(x), label(class_label), weight(w)
  {}


  KolmogorovSmirnov::KolmogorovSmirnov(void)
    : cached_(true), score_(0.0), sum_w1_(0.0), sum_w2_(0.0)
  {
  }


  void KolmogorovSmirnov::add(double x, bool target, double w)
  {
    if (w==0) // ignoring zero weight data
      return;

    assert(!std::isnan(x));
    data_.insert(Element(x,target,w));
    if (target){
      sum_w1_+=w;
    }
    else {
      sum_w2_+=w;
    }
    cached_=false;
  }


  double KolmogorovSmirnov::p_value(void) const
  {
    double res=0;
    double res2=0;
    double s2 = score() * score() * sum_w1_*sum_w2_/(sum_w1_+sum_w2_);
    int sign = 1;
    for (size_t k = 1; k<100; ++k) {
      res += sign * 2 * exp(-2.0 * k * k * s2);
      sign *= -1;
      // ignore remaining terms as they are small
      if (res==res2)
        return res;
      res2 = res;
    }

    return res;
  }


  double KolmogorovSmirnov::p_value(size_t perm) const
  {
    size_t count=0;
    double score_threshold = score()-10*std::numeric_limits<double>().epsilon();
    KolmogorovSmirnov ks(*this);

    for (size_t i=0; i<perm; ++i){
      ks.shuffle();
      if (ks.score()>=score_threshold)
        ++count;
    }
    return static_cast<double>(count)/perm;
  }


  void KolmogorovSmirnov::remove(double value, bool class_label, double weight)
  {
    if (!weight)
      return;
    Element e(value, class_label, weight);
    typedef std::multiset<Element>::iterator iterator;
    std::pair<iterator, iterator> iter = data_.equal_range(e);
    while (iter.first!=iter.second) {
      if (iter.first->label==class_label && iter.first->weight==weight) {
        if (class_label)
          sum_w1_-=weight;
        else
          sum_w2_-=weight;
        data_.erase(iter.first);
        cached_=false;
        return;
      }
      ++iter.first;
    }
    std::ostringstream ss;
    ss << "KolmogorovSmirnov::remove: " << value << " " << class_label
       << " " << weight;
    throw utility::runtime_error(ss.str());
  }


  void KolmogorovSmirnov::reset(void)
  {
    cached_=false;
    data_.clear();
    sum_w1_=0;
    sum_w2_=0;
  }


  double KolmogorovSmirnov::score(void) const
  {
    return std::abs(signed_score());
  }


  void KolmogorovSmirnov::scores(std::vector<double>& res) const
  {
    res.clear();
    res.reserve(data_.size());
    data_w::const_iterator iter(data_.begin());
    double f1=0;
    double f2=0;
    while(iter!=data_.end()){
      size_t count=0;
      double value = iter->value;
      while (iter!=data_.end() && iter->value==value) {
        if (iter->label) {
          f1 += iter->weight;
        }
        else {
          f2 += iter->weight;
        }
        ++count;
        ++iter;
      }
      res.resize(res.size()+count, f1/sum_w1_-f2/sum_w2_);
    }
  }


  void KolmogorovSmirnov::shuffle(void)
  {
    std::deque<bool> labels;
    for (data_w::const_iterator iter=data_.begin();
         iter!=data_.end(); ++iter) {
      labels.push_back(iter->label);
    }
    random::random_shuffle(labels.begin(), labels.end());
    std::deque<bool>::const_iterator label = labels.begin();
    for (data_w::iterator iter=data_.begin();
         iter!=data_.end(); ++iter, ++label) {
      // label does not affect sorting of set, so modifying label is safe
#ifdef MULTISET_ITERATOR_IS_MUTABLE
      (*iter).label = *label;
#else
      // standard requires that multiset::iterator is modifiable, but
      // some implementations do not fulfill the requirement
      const_cast<Element&>(*iter).label = *label;
#endif
    }
    sum_w1_ = 0;
    sum_w2_ = 0;
    add_sum_w(data_.begin(), data_.end());
    cached_ = false;
  }


  double KolmogorovSmirnov::signed_score(void) const
  {
    if (cached_)
      return score_;
    if (!sum_w1_ || !sum_w2_)
      return 0.0;
    score_=0;
    std::vector<double> v;
    scores(v);
    std::less<double> l;
    utility::abs<double> a;
    using utility::make_compose_f_gx_hy;
    score_ = *std::max_element(v.begin(), v.end(),
                               make_compose_f_gx_hy(l, a, a));
    cached_=true;
    return score_;
  }


  std::ostream& operator<<(std::ostream& os, const KolmogorovSmirnov& ks)
  {
    std::vector<double> s;
    ks.scores(s);
    for (size_t i=0; i<s.size(); ++i)
      os << i << "\t" << s[i] << "\n";
    return os;
  }


  bool KolmogorovSmirnov::Element::operator<
  (const KolmogorovSmirnov::Element& rhs) const
  {
    return value<rhs.value;
  }

}}} // of namespace theplu yat statistics