source: trunk/yat/statistics/KolmogorovSmirnovOneSample.cc @ 2998

// $Id: KolmogorovSmirnovOneSample.cc 2998 2013-03-14 08:08:54Z peter $

/*
  Copyright (C) 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 "KolmogorovSmirnovOneSample.h"

#include "yat/utility/Exception.h"

#include <cassert>
#include <cmath>
#include <sstream>

#include <iostream>

namespace theplu {
namespace yat {
namespace statistics {

  KolmogorovSmirnovOneSample::KolmogorovSmirnovOneSample(void)
    : sum_w_(0)
  {
  }


  void KolmogorovSmirnovOneSample::add(double value, double weight)
  {
    if (weight==0)
      return;
    cached_ = false;
    data_.insert(std::make_pair(value, weight));
    sum_w_ += weight;
  }


  double KolmogorovSmirnovOneSample::p_value(void) const
  {
    double res=0;
    double res2=0;
    double s2 = score() * score() * sum_w_;
    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;
  }


  void KolmogorovSmirnovOneSample::remove(double value, double weight)
  {
    if (!weight)
      return;
    typedef std::multimap<double, double>::iterator iterator;
    std::pair<iterator, iterator> range = data_.equal_range(value);
    for ( ; range.first!=range.second; ++range.first) {
      if (range.first->second == weight) {
        data_.erase(range.first);
        sum_w_ -= weight;
        cached_ = false;
        return;
      }
    }
    std::ostringstream ss;
    ss << "KolmogorovSmirnovOneSample::remove: " << value << " " << weight;
    throw utility::runtime_error(ss.str());
  }


  void KolmogorovSmirnovOneSample::reset(void)
  {
    data_.clear();
    sum_w_ = 0;
  }


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


  double KolmogorovSmirnovOneSample::signed_score(void) const
  {
    if (cached_)
      return score_;
    score_ = 0.0;
    if (data_.empty())
      return 0.0;

    typedef std::multimap<double, double>::const_iterator iterator;
    iterator iter = data_.begin();
    double w = 0;
    while (iter != data_.end()) {
      double x = iter->first;
      while (iter!=data_.end() && iter->first==x) {
        w += iter->second;
        ++iter;
      }
      double s = w / sum_w_ - x;
      if (std::abs(s) > std::abs(score_))
        score_ = s;
    }

    cached_ = true;
    return score_;
  }

}}} // of namespace theplu yat statistics