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

#ifndef _theplu_yat_statistics_utility_ 
#define _theplu_yat_statistics_utility_ 

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

/*
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005 Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2007, 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 "Percentiler.h"

#include "yat/classifier/DataLookupWeighted1D.h"
#include "yat/classifier/Target.h"
#include "yat/utility/deprecate.h"
#include "yat/utility/iterator_traits.h"
#include "yat/utility/VectorBase.h"
#include "yat/utility/yat_assert.h"

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

#include <gsl/gsl_statistics_double.h>

namespace theplu {
namespace yat {
namespace statistics {  

  /**
     Adding a range [\a first, \a last) into an object of type T. The
     requirements for the type T is to have an add(double, bool, double)
     function.
  */
  template <typename T, typename ForwardIterator>
  void add(T& o, ForwardIterator first, ForwardIterator last,
           const classifier::Target& target);

  ///
  /// 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(unsigned int k, unsigned int n1, 
                              unsigned int n2, unsigned int t);


  /**
     \brief one-sided p-value

     This function uses the t-distribution to calculate the one-sided
     p-value. Given that the true correlation is zero (Null
     hypothesis) the estimated correlation, r, after a transformation
     is t-distributed:

     \f$ \sqrt{(n-2)} \frac{r}{\sqrt{(1-r^2)}} \in t(n-2) \f$

     \return Probability that correlation is larger than \a r by
     chance when having \a n samples. For \a r larger or equal to 1.0,
     0.0 is returned. For \a r smaller or equal to -1.0, 1.0 is
     returned.
   */
  double pearson_p_value(double r, unsigned int n);

  ///
  /// @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::VectorBase&);


  ///
  /// @brief Median absolute deviation from median
  ///
  /// Function is non-mutable function
  ///
  template <class T>
  double mad(T first, T last, const bool sorted=false);


  ///
  /// 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); 


  /**
     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: RandomAccessIterator is an iterator over a range of
     doubles (or any type being convertable to double).
     
     @return \a p'th percentile of range

     \deprecated percentile2 will replace this function in the future

     \note the definition of percentile used here is not identical to
     that one used in percentile2 and Percentile. The difference is
     smaller for large ranges.
  */
  template <class RandomAccessIterator>
  double percentile(RandomAccessIterator first, RandomAccessIterator last, 
                    double p, bool sorted=false) YAT_DEPRECATE;
  

  /**
     \see Percentiler
     
     \since new in yat 0.5
   */
  template <class RandomAccessIterator>
  double percentile2(RandomAccessIterator first, RandomAccessIterator last, 
                     double p, bool sorted=false)
  {
    Percentiler percentiler(p, sorted);
    return percentiler(first, last);
  }

  ///
  /// @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::VectorBase&);
  

  //////// template implementations ///////////

  template <typename T, typename ForwardIterator>
  void add(T& o, ForwardIterator first, ForwardIterator last,
           const classifier::Target& target)
  {
    for (size_t i=0; first!=last; ++i, ++first)
      o.add(utility::iterator_traits<ForwardIterator>().data(first),
            target.binary(i), 
            utility::iterator_traits<ForwardIterator>().weight(first));
  } 


  template <class T>
  double mad(T first, T last, const bool sorted=false)
  {
    double m = median(first, last, sorted);
    typedef std::vector<typename std::iterator_traits<T>::value_type> Vec;
    Vec ad;
    ad.reserve(std::distance(first, last));
    // copy values (and weights if T is weighted)
    std::copy(first, last, std::back_insert_iterator<Vec>(ad));
    utility::iterator_traits<typename Vec::iterator> traits;
    for (typename Vec::iterator i=ad.begin() ; i!=ad.end(); ++i) {
      // adjust data
      traits.data(i) = std::abs(traits.data(i)-m);
    }
    std::sort(ad.begin(), ad.end());
    return median(ad.begin(), ad.end(), true);
  }
  

  template <class T> 
  double median(T first, T last, const bool sorted=false) 
  { 
    return percentile2(first, last, 50.0, sorted); 
  }


  template <class RandomAccessIterator>
  double percentile(RandomAccessIterator first, RandomAccessIterator last, 
                    double p, bool sorted=false)
  {
    // range is one value only is a special case
    if (first+1 == last)
      return utility::iterator_traits<RandomAccessIterator>().data(first);
    if (sorted) {
      // have to take care of this special case
      if (p>=100)
        return utility::iterator_traits<RandomAccessIterator>().data(--last);
      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<typename std::iterator_traits<RandomAccessIterator>::value_type> 
      v_copy;
    v_copy.reserve(std::distance(first,last));
    std::copy(first, last, std::back_inserter(v_copy));
    size_t i = static_cast<size_t>(p/100 * (v_copy.size()-1));
    if (i+2 < v_copy.size()) {
      std::partial_sort(v_copy.begin(), v_copy.begin()+i+2, v_copy.end());
    }
    else
      std::sort(v_copy.begin(), v_copy.end());
    return percentile(v_copy.begin(), v_copy.end(), p, true);
  }


}}} // of namespace statistics, yat, and theplu

#endif