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

// $Id: KolmogorovSmirnov.cc 1437 2008-08-25 17:55:00Z peter $

/*
  Copyright (C) 2008 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 2 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 this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

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

#include <algorithm>
#include <cassert>
#include <cmath>
#include <deque>
#include <functional>


namespace theplu {
namespace yat {
namespace statistics {


  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;

    data_.insert(std::make_pair(std::make_pair(x,w),target));
    if (target){
      ++n1_;
      sum_w1_+=w;
    }
    else {
      ++n2_;
      sum_w2_+=w;
    }
    cached_=false;
  }


  double KolmogorovSmirnov::p_value(size_t perm) const
  {
    size_t count=0;
    std::deque<bool> targets;
    typedef data_w::const_iterator const_iter;
    for (const_iter i(data_.begin()); i!=data_.end(); ++i)
      targets.push_back(i->second);

    for (size_t i=0; i<perm; ++i){
      random::random_shuffle(targets.begin(), targets.end()); 
      KolmogorovSmirnov ks;
      std::deque<bool>::const_iterator target_i(targets.begin());
      const_iter end(data_.end());
      for (const_iter i(data_.begin()); i!=end; ++i, ++target_i) 
        ks.add(i->first.first, *target_i, i->first.second);
      
      if (ks.score()>=score())
        ++count;
    }
    return static_cast<double>(count)/perm;
  }


  double KolmogorovSmirnov::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::abs(*std::max_element(v.begin(), v.end(), 
                                        make_compose_f_gx_hy(l, a, a)));
    return 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;
    // Peter, we should take care of ties!
    while(iter!=data_.end()){
      if (iter->second)
        f1 += iter->first.second;
      else
        f2 += iter->first.second;
      res.push_back(f1/sum_w1_-f2/sum_w2_);
      ++iter;
    }
  }


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


  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;
  }


}}} // of namespace theplu yat statistics