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

#ifndef _theplu_yat_statistics_utility_
#define _theplu_yat_statistics_utility_

// $Id: utility.h 3245 2014-05-24 17:01:31Z 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, 2010, 2011, 2013, 2014 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/normalizer/utility.h"
#include "yat/utility/concept_check.h"
#include "yat/utility/DataIterator.h"
#include "yat/utility/deprecate.h"
#include "yat/utility/iterator_traits.h"
#include "yat/utility/sort_index.h"
#include "yat/utility/Vector.h"
#include "yat/utility/VectorBase.h"
#include "yat/utility/yat_assert.h"

#include <boost/concept_check.hpp>
#include <boost/iterator/permutation_iterator.hpp>

#include <gsl/gsl_math.h>
#include <gsl/gsl_statistics_double.h>

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

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

  /**
     \brief Benjamini Hochberg multiple test correction

     Given a sorted range of p-values such that
     \f$ p_1 \le p_2 \le ... \le p_N \f$ a Benjamnini-Hochberg corrected
     p-value, \c q, is calculated recursively as
     \f$ q_i = \f$ min \f$(p_i \frac{N}{i}, q_{i+1})\f$ with the anchor
     constraint that \f$ q_m = p_m \f$.

     Type Requirements:
     - \c BidirectionalIterator1 is a \bidirectional_iterator
     - \c BidirectionalIterator1::value type is convertible to \c double
     - \c BidirectionalIterator2 is a mutable \bidirectional_iterator
     - \c double is convertible to \c BidirectionalIterator2::value_type

     \since New in yat 0.8
   */
  template<typename BidirectionalIterator1, typename BidirectionalIterator2>
  void benjamini_hochberg(BidirectionalIterator1 first,
                          BidirectionalIterator1 last,
                          BidirectionalIterator2 result);


  /**
     \brief Benjamini Hochberg multiple test correction

     Similar to benjamini_hochberg() but does not assume that input
     range, [first, last), is sorted. The resulting range is the same
     as if sorting input range, call benjamini_hochberg, and unsort
     the result range.

     Type Requirements:
     - \c RandomAccessIterator models a \random_access_iterator
     - \c RandomAccessIterator::value type is convertible to \c double
     - \c MutableRandomAccessIterator models a mutable \random_access_iterator
     - \c double is convertible to \c MutableRandomAccessIterator::value_type

     \since New in yat 0.13
   */
  template<typename RandomAccessIterator, typename MutableRandomAccessIterator>
  void benjamini_hochberg_unsorted(RandomAccessIterator first,
                                   RandomAccessIterator last,
                                   MutableRandomAccessIterator result);


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


  /**
     The entropy is calculated as \f$ - \sum_i p_i \log p_i \f$ where
     \f$p_i = \frac{n_i}{\sum_j n_j} \f$

     Requirements:
     - \c InputIterator should be an \input_iterator
     - \c InputIterator::value_type must be convertible to \c double

     \since New in yat 0.12
   */
  template<typename InputIterator>
  double entropy(InputIterator first, InputIterator last);

  /**
     \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
  ///
  /// Requirements: \c RandomAccessIterator should be a \ref
  /// concept_data_iterator and \random_access_iterator
  ///
  /// Since 0.6 function also work with a \ref
  /// concept_weighted_iterator
  ///
  template <class RandomAccessIterator>
  double mad(RandomAccessIterator first, RandomAccessIterator last,
             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: \c RandomAccessIterator should be a \ref
  /// concept_data_iterator and \random_access_iterator
  ///
  /// @return median of range
  ///
  template <class RandomAccessIterator>
  double median(RandomAccessIterator first, RandomAccessIterator last,
                bool sorted=false);


  /**
     \brief Calculates the mutual information of \a A.

     The elements in A are unnormalized probabilies of the joint
     distribution.

     The mutual information is calculated as \f$ \sum \sum p(x,y) \log_2
     \frac {p(x,y)} {p(x)p(y)} \f$ where
     \f$ p(x,y) = \frac {A_{xy}}{\sum_{x,y} A_{xy}} \f$;
     \f$ p(x) = \sum_y A_{xy} / \sum_{x,y} A_{xy} \f$;
     \f$ p(y) = \sum_x A_{xy} / \sum_{x,y} A_{xy} \f$

     Requirements:
     - \c T must be a model of \ref concept_container_2d
     - \c T::value_type must be convertible to \c double

     \return mutual information in bits; if you want in natural base
     multiply with \c M_LN2 (defined in \c gsl/gsl_math.h )

     \since New in yat 0.12
   */
  template<class T>
  double mutual_information(const T& A);

  /**
     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)
  {
    utility::iterator_traits<ForwardIterator> traits;
    for (size_t i=0; first!=last; ++i, ++first)
      o.add(traits.data(first), target.binary(i), traits.weight(first));
  }


  template<typename BidirectionalIterator1, typename BidirectionalIterator2>
  void benjamini_hochberg(BidirectionalIterator1 first,
                          BidirectionalIterator1 last,
                          BidirectionalIterator2 result)
  {
    using boost::Mutable_BidirectionalIterator;
    using boost::BidirectionalIterator;
    BOOST_CONCEPT_ASSERT((BidirectionalIterator<BidirectionalIterator1>));
    BOOST_CONCEPT_ASSERT((Mutable_BidirectionalIterator<BidirectionalIterator2>));
    typedef typename std::iterator_traits<BidirectionalIterator1> traits;
    typename traits::difference_type n = std::distance(first, last);
    if (!n)
      return;
    std::advance(result, n-1);
    first = last;
    --first;
    typename traits::difference_type rank = n;

    double prev = 1.0;
    while (rank) {
      *result = std::min(*first * n/static_cast<double>(rank), prev);
      prev = *result;
      --rank;
      --first;
      --result;
    }
  }


  template<typename RandomAccessIterator, typename MutableRandomAccessIterator>
  void benjamini_hochberg_unsorted(RandomAccessIterator first,
                                   RandomAccessIterator last,
                                   MutableRandomAccessIterator result)
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
    using boost::Mutable_RandomAccessIterator;
    BOOST_CONCEPT_ASSERT((Mutable_RandomAccessIterator<MutableRandomAccessIterator>));

    std::vector<size_t> idx;
    utility::sort_index(first, last, idx);
    benjamini_hochberg(boost::make_permutation_iterator(first, idx.begin()),
                       boost::make_permutation_iterator(first, idx.end()),
                       boost::make_permutation_iterator(result, idx.begin()));
  }


  template<typename InputIterator>
  double entropy(InputIterator first, InputIterator last)
  {
    BOOST_CONCEPT_ASSERT((boost::InputIterator<InputIterator>));
    using boost::Convertible;
    typedef typename InputIterator::value_type T;
    BOOST_CONCEPT_ASSERT((Convertible<T,double>));
    double sum = 0;
    double N = 0;
    for (; first != last; ++first) {
      if (*first) {
        N += *first;
        sum += *first * std::log(static_cast<double>(*first));
      }
    }
    return -sum / N + log(N);
  }


  template <class RandomAccessIterator>
  double mad(RandomAccessIterator first, RandomAccessIterator last,
             bool sorted)
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
    double m = median(first, last, sorted);
    typedef typename std::iterator_traits<RandomAccessIterator>::value_type T;
    std::vector<T> ad(std::distance(first, last));
    // copy weights if needed
    normalizer::detail::copy_weight_if_weighted(first, last, ad.begin());
    // assign data (absolute deviation from median)
    utility::iterator_traits<RandomAccessIterator> traits;
    utility::DataIterator<typename std::vector<T>::iterator> first2(ad.begin());
    while (first!=last) {
      *first2 = std::abs(traits.data(first)-m);
      ++first;
      ++first2;
    }
    return median(ad.begin(), ad.end(), false);
  }


  template <class RandomAccessIterator>
  double median(RandomAccessIterator first, RandomAccessIterator last,
                bool sorted)
  {
    return percentile2(first, last, 50.0, sorted);
  }


  template<class T>
  double mutual_information(const T& n)
  {
    BOOST_CONCEPT_ASSERT((utility::Container2D<T>));
    using boost::Convertible;
    BOOST_CONCEPT_ASSERT((Convertible<typename T::value_type,double>));

    // p_x = \sum_y p_xy

    // Mutual Information is defined as
    // \sum_xy p_xy * log (p_xy / (p_x p_y)) =
    // \sum_xy p_xy * [log p_xy - log p_x - log p_y]
    // \sum_xy p_xy log p_xy - p_xy log p_x - p_xy log p_y
    // \sum_xy p_xy log p_xy - \sum_x p_x log p_x - \sum_y p_y log p_y
    // - entropy_xy + entropy_x + entropy_y

    utility::Vector rowsum(n.columns(), 0);
    for (size_t c = 0; c<n.columns(); ++c)
      rowsum(c) = std::accumulate(n.begin_column(c), n.end_column(c), 0);

    utility::Vector colsum(n.rows(), 0);
    for (size_t r = 0; r<n.rows(); ++r)
      colsum(r) = std::accumulate(n.begin_row(r), n.end_row(r), 0);

    double mi  = - entropy(n.begin(), n.end());
    mi += entropy(rowsum.begin(), rowsum.end());
    mi += entropy(colsum.begin(), colsum.end());

    return mi/M_LN2;
  }


  template <class RandomAccessIterator>
  double percentile(RandomAccessIterator first, RandomAccessIterator last,
                    double p, bool sorted)
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
    // 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];
    }
    using std::iterator_traits;
    typedef typename iterator_traits<RandomAccessIterator>::value_type value_t;
    std::vector<value_t> v_copy(first, last);
    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