utility.h

#ifndef _theplu_yat_statistics_utility_ 
#define _theplu_yat_statistics_utility_ 

// $Id: utility.h 2625 2011-11-08 00:07:24Z 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 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/VectorBase.h"
#include "yat/utility/yat_assert.h"

#include <boost/concept_check.hpp>

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

#include <gsl/gsl_statistics_double.h>

namespace theplu {
namespace yat {
namespace statistics {  

  template <typename T, typename ForwardIterator>
  void add(T& o, ForwardIterator first, ForwardIterator last,
           const classifier::Target& target);

  template<typename BidirectionalIterator1, typename BidirectionalIterator2>
  void benjamini_hochberg(BidirectionalIterator1 first,
                          BidirectionalIterator1 last,
                          BidirectionalIterator2 result);
                          

  double cdf_hypergeometric_P(unsigned int k, unsigned int n1, 
                              unsigned int n2, unsigned int t);


  double pearson_p_value(double r, unsigned int n);

  double kurtosis(const utility::VectorBase&);


  template <class RandomAccessIterator>
  double mad(RandomAccessIterator first, RandomAccessIterator last, 
             bool sorted=false);


  template <class RandomAccessIterator> 
  double median(RandomAccessIterator first, RandomAccessIterator last, 
                bool sorted=false); 


  template <class RandomAccessIterator>
  double percentile(RandomAccessIterator first, RandomAccessIterator last, 
                    double p, bool sorted=false) YAT_DEPRECATE;

  template <class RandomAccessIterator>
  double percentile2(RandomAccessIterator first, RandomAccessIterator last, 
                     double p, bool sorted=false)
  {
    Percentiler percentiler(p, sorted);
    return percentiler(first, last);
  }

  double skewness(const utility::VectorBase&);
  


  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;
    BOOST_CONCEPT_ASSERT((boost::BidirectionalIterator<BidirectionalIterator1>));
    BOOST_CONCEPT_ASSERT((Mutable_BidirectionalIterator<BidirectionalIterator2>));
    size_t n = std::distance(first, last);
    if (!n)
      return;
    std::advance(result, n-1);
    first = last;
    --first;
    size_t rank = n;
    
    double prev = 1.0;
    while (rank) {
      *result = std::min(*first * n/static_cast<double>(rank), prev);
      prev = *result;
      --rank;
      --first;
      --result;
    }
  }


  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;
    }
    std::sort(ad.begin(), ad.end());
    return median(ad.begin(), ad.end(), true);
  }
  

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


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