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

#ifndef _theplu_yat_statistics_kolmogorov_smirnov_ 
#define _theplu_yat_statistics_kolmogorov_smirnov_ 

// $Id: KolmogorovSmirnov.h 1874 2009-03-17 22:09:01Z peter $

/*
  Copyright (C) 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009 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 <set>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {  

  /**
     \brief Kolmogow 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
      */
      inline Element(void) {};

      /**
         \brief Constructor
       */
      inline Element(double x, bool class_label, double w=1.0)
        : value(x), label(class_label), weight(w) {};

      /**
         \brief data value
       */
      double value;
      
      /** 
          bool telling which group the data point belong 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 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
   */
  std::ostream& operator<<(std::ostream&, const KolmogorovSmirnov&);


  //  template implementations

  template <typename ForwardIterator>
  void KolmogorovSmirnov::add(ForwardIterator first, ForwardIterator last)
  {
    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