source: trunk/yat/statistics/utility.h @ 934

#ifndef _theplu_yat_statistics_utility_ 
#define _theplu_yat_statistics_utility_ 

// $Id: utility.h 934 2007-10-05 21:27:34Z peter $

/*
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005 Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
  Copyright (C) 2007 Jari Häkkinen, Peter Johansson

  This file is part of the yat library, http://trac.thep.lu.se/trac/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 "yat/classifier/DataLookupWeighted1D.h"
#include "yat/classifier/Target.h"
#include "yat/utility/vector.h"
#include "yat/utility/yat_assert.h"

#include <algorithm>
#include <cmath>
#include <vector>

#include <gsl/gsl_statistics_double.h>

namespace theplu {
namespace yat {
namespace statistics {  

  //forward declarations
  template <class T> 
  double median(T first, T last, const bool sorted=false); 

  template <class T>
  double percentile(T first, T last, double p, bool sorted=false);
  
  /**
     Adding each value in an array \a v \a to an object \a o.  The
     requirements for the type T1 is to have an add(double, bool)
     function, and for T2 of the array \a v are: operator[] returning
     an element and function size() returning the number of elements.
  */
  template <typename T1, typename T2>
  void add(T1& o, const T2& v, const classifier::Target& target)
  {
    for (size_t i=0; i<v.size(); ++i) 
      o.add(v[i],target.binary(i));
  } 

  /**
     Adding each value in an array \a v \a to an object \a o.  The
     requirements for the type T1 is to have an add(double, bool)
     function, and for T2 of the array \a v are: operator[] returning
     an element and function size() returning the number of elements.
  */
  template <typename T1>
  void add(T1& o, const classifier::DataLookupWeighted1D& v, 
           const classifier::Target& target)
  {
    for (size_t i=0; i<v.size(); ++i) 
      o.add(v.data(i),target.binary(i),v.weight(i));
  } 

  ///
  /// Calculates the probability to get \a k or smaller from a
  /// hypergeometric distribution with parameters \a n1 \a n2 \a
  /// t. Hypergeomtric situation you get in the following situation:
  /// Let there be \a n1 ways for a "good" selection and \a n2 ways
  /// for a "bad" selection out of a total of possibilities. Take \a
  /// t samples without replacement and \a k of those are "good"
  /// samples. \a k will follow a hypergeomtric distribution.
  ///
  /// @return cumulative hypergeomtric distribution functions P(k).
  ///
  double cdf_hypergeometric_P(u_int k, u_int n1, u_int n2, u_int t);


  ///
  /// @brief Computes the kurtosis of the data in a vector.
  ///
  /// The kurtosis measures how sharply peaked a distribution is,
  /// relative to its width. The kurtosis is normalized to zero for a
  /// gaussian distribution.
  ///
  double kurtosis(const utility::vector&);


  ///
  /// @brief Median absolute deviation from median
  ///
  /// Function is non-mutable function
  ///
  template <class T>
  double mad(T first, T last, const bool sorted=false)
  {
    double m = median(first, last, sorted);
    std::vector<double> ad;
    ad.reserve(std::distance(first, last));
    for( ; first!=last; ++first)
      ad.push_back(fabs(*first-m));
    std::sort(ad.begin(), ad.end());
    return median(ad.begin(), ad.end(), true);
  }
  

  ///
  /// Median is defined to be value in the middle. If number of values
  /// is even median is the average of the two middle values.  the
  /// median value is given by p equal to 50. If \a sorted is false
  /// (default), the range is copied, the copy is sorted, and then
  /// used to calculate the median.
  ///
  /// Function is a non-mutable function, i.e., \a first and \a last
  /// can be const_iterators.
  ///
  /// Requirements: T should be an iterator over a range of doubles (or
  /// any type being convertable to double). 
  ///
  /// @return median of range
  ///
  template <class T> 
  double median(T first, T last, const bool sorted=false) 
  { return percentile(first, last, 50.0, sorted); }

  /**
     The percentile is determined by the \a p, a number between 0 and
     100. The percentile is found by interpolation, using the formula
     \f$ percentile = (1 - \delta) x_i + \delta x_{i+1} \f$ where \a
     p is floor\f$((n - 1)p/100)\f$ and \f$ \delta \f$ is \f$
     (n-1)p/100 - i \f$.Thus the minimum value of the vector is given
     by p equal to zero, the maximum is given by p equal to 100 and
     the median value is given by p equal to 50. If @a sorted
     is false (default), the vector is copied, the copy is sorted,
     and then used to calculate the median.

     Function is a non-mutable function, i.e., \a first and \a last
     can be const_iterators.
     
     Requirements: T should be an iterator over a range of doubles (or
     any type being convertable to double). If \a sorted is false
     iterator must be mutable, else read-only iterator is also ok.
     
     @return \a p'th percentile of range
  */
  template <class T>
  double percentile(T first, T last, double p, bool sorted=false)
  {
    yat_assert(first<last && "range is invalid");
    yat_assert(p>=0);
    yat_assert(p<=100);
    if (sorted){
      if (p>=100)
        return *(--last);
      // range is one value only is a special case
      if (std::distance(first, last)==1) 
        return *first;
      double j = p/100 * (std::distance(first,last)-1);
      int i = static_cast<int>(j);
      return (1-j+floor(j))*first[i] + (j-floor(j))*first[i+1];
    }

    std::vector<double> v_copy;
    v_copy.reserve(std::distance(first,last));
    std::copy(first, last, std::back_inserter(v_copy));
    double j = p/100 * (v_copy.size()-1);
    int i = static_cast<int>(j);
    std::partial_sort(v_copy.begin(),v_copy.begin()+i+2 , v_copy.end());
    return percentile(v_copy.begin(), v_copy.end(), p, true);
  }

  ///
  /// @brief Computes the skewness of the data in a vector.
  ///
  /// The skewness measures the asymmetry of the tails of a
  /// distribution.
  ///
  double skewness(const utility::vector&);
  
}}} // of namespace statistics, yat, and theplu

#endif