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

#ifndef _theplu_yat_statistics_kolmogorov_smirnov_
#define _theplu_yat_statistics_kolmogorov_smirnov_

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

/*
  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2010, 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 <boost/concept_check.hpp>

#include <iosfwd>
#include <set>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {

  /**
     \brief Kolmogorov Smirnov Test
   */
  class KolmogorovSmirnov
  {
  public:
    /**
       struct used to store data in KolmogorovSmirnov.

       This struct is public to make usage of the add range function
       more convenient.

       \since New in yat 0.5
     */
    struct Element
    {
      /**
         \brief default constructor
      */
      Element(void);

      /**
         \brief Constructor
       */
      Element(double x, bool class_label, double w=1.0);

      /**
         \brief data value
       */
      double value;

      /**
          bool telling which group the data point belongs to
      */
      bool label;

      /**
         weight for the data point
       */
      double weight;

      /**
         \return true if value is less than rhs.value
       */
      bool operator<(const Element& rhs) const;
    };

    /**
       \brief Constructor
     */
    KolmogorovSmirnov(void);

    /**
       \brief add a value
     */
    void add(double value, bool class_label, double weight=1.0);

    /**
       \brief add a range

       value_type of \a ForwardIterator must be convertible to
       KolmogorovSmirnov::Element

       Insertion takes typically N*log(N). However, if range is
       sorted, insertion takes linear time. A common use case of this
       function is when calculating several KS statistics on a data
       set (values and weights) with different class labels.

       \since New in yat 0.5
    */
    template <typename ForwardIterator>
    void add(ForwardIterator first, ForwardIterator last);

    /**
       \brief Large-Sample Approximation

       This analytical approximation of p-value can be used when all
       weight equal unity and sample sizes \a n and \a m are
       large. The p-value is calcuated as \f$ P = \displaystyle - 2
       \sum_{k=1}^{\infty} (-1)^ke^{-2k^2s^2}\f$, where s is the
       scaled score:

       \f$ s = \sqrt\frac{nm}{n+m} \f$ score().

       \since New in yat 0.5

       Following Hollander and Wolfe
    */
    double p_value(void) const;

    /**
       \brief p-value

       Performs a permutation test using \a perm label randomizations
       and calculate how often a score equal or larger than score() is
       obtained.

       \see shuffle
    */
    double p_value(size_t perm) const;

    /**
       \brief Remove a data point

       \throw utility::runtime_error if no data point exist with \a
       value, \a class_label, and \a weight.

       \since New in yat 0.9
     */
    void remove(double value, bool class_label, double weight=1.0);

    /**
       \brief resets everything to zero
    */
    void reset(void);

    /**
       \brief Kolmogorov Smirnov statistic

       \f$ sup_x | F_{\textrm{True}}(x) - F_{\textrm{False}}(x) | \f$ where
       \f$ F(x) = \frac{\sum_{i:x_i\le x}w_i}{ \sum w_i} \f$
    */
    double score(void) const;

    /**
       \brief shuffle class labels

       This is equivalent to reset and re-add values with shuffled
       class labels.

       \since New in yat 0.5
     */
    void shuffle(void);

    /**
       Same as score() but keeping the sign, in other words,
       abs(signed_score())==score()

       A positive score implies that values in class \c true \em on
       \em average are smaller than values in class \c false.

       \since New in yat 0.5
     */
    double signed_score(void) const;

  private:
    void scores(std::vector<double>&) const;
    // add weights to sum_w1 and sum_w2 respectively depending on
    // label in element.
    template <typename ForwardIterator>
    void add_sum_w(ForwardIterator first, ForwardIterator last);

    mutable bool cached_;
    mutable double score_;
    typedef std::multiset<Element> data_w;
    data_w data_;
    double sum_w1_;
    double sum_w2_;

    friend std::ostream& operator<<(std::ostream&, const KolmogorovSmirnov&);

    // using compiler generated copy and assignment
    //KolmogorovSmirnov(const KolmogorovSmirnov&);
    //KolmogorovSmirnov& operator=(const KolmogorovSmirnov&);
  };

  /**
     \brief output operator

     \relates KolmogorovSmirnov
   */
  std::ostream& operator<<(std::ostream&, const KolmogorovSmirnov&);


  //  template implementations

  template <typename ForwardIterator>
  void KolmogorovSmirnov::add(ForwardIterator first, ForwardIterator last)
  {
    BOOST_CONCEPT_ASSERT((boost::ForwardIterator<ForwardIterator>));
    typedef typename std::iterator_traits<ForwardIterator>::reference ref;
    BOOST_CONCEPT_ASSERT((boost::Convertible<ref, KolmogorovSmirnov::Element>));
    ForwardIterator iter(first);
    typename data_w::const_iterator hint(data_.begin());
    for ( ; iter!=last; ++iter)
      if ((*iter).weight) // ignore data points with zero weight
        hint = data_.insert(hint, *iter);
    add_sum_w(first, last);
    cached_=false;
  }


  template <typename ForwardIterator>
  void KolmogorovSmirnov::add_sum_w(ForwardIterator first,
                                    ForwardIterator last)
  {
    while (first!=last) {
      if ((*first).label)
        sum_w1_ += (*first).weight;
      else
        sum_w2_ += (*first).weight;
      ++first;
    }
  }

}}} // of namespace theplu yat statistics

#endif