Changeset 3236 for trunk/yat/statistics/utility.h

Timestamp:

May 23, 2014, 3:42:51 PM (9 years ago)

Author:

Peter

Message:

BH correction on unsorted range. closes #796

File:

• trunk/yat/statistics/utility.h (modified) (12 diffs)

trunk/yat/statistics/utility.h

@@ -9 +9 @@
   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 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
@@ -36 +36 @@
 #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"
@@ -41 +42 @@
 
 #include <boost/concept_check.hpp>
+#include <boost/iterator/permutation_iterator.hpp>
+
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_statistics_double.h>
@@ -72 +75 @@
      constraint that \f$ q_m = p_m \f$.
 
-     Requirements: \c BidirectionalIterator1 should be a
-     \bidirectional_iterator and \c BidirectionalIterator2 should be a
-     mutable \bidirectional_iterator
+     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
@@ -82 +87 @@
                           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);
 
 
@@ -95 +122 @@
   /// @return cumulative hypergeomtric distribution functions P(k).
   ///
-  double cdf_hypergeometric_P(unsigned int k, unsigned int n1, 
+  double cdf_hypergeometric_P(unsigned int k, unsigned int n1,
                               unsigned int n2, unsigned int t);
 
@@ -248 +275 @@
   ///
   double skewness(const utility::VectorBase&);
-  
+
 
   //////// template implementations ///////////
@@ -259 +286 @@
     for (size_t i=0; first!=last; ++i, ++first)
       o.add(traits.data(first), target.binary(i), traits.weight(first));
-  } 
+  }
 
 
@@ -268 +295 @@
   {
     using boost::Mutable_BidirectionalIterator;
-    BOOST_CONCEPT_ASSERT((boost::BidirectionalIterator<BidirectionalIterator1>));
+    using boost::BidirectionalIterator;
+    BOOST_CONCEPT_ASSERT((BidirectionalIterator<BidirectionalIterator1>));
     BOOST_CONCEPT_ASSERT((Mutable_BidirectionalIterator<BidirectionalIterator2>));
-    size_t n = std::distance(first, last);
+    typedef typename std::iterator_traits<BidirectionalIterator1> traits;
+    typename traits::difference_type n = std::distance(first, last);
     if (!n)
       return;
@@ -276 +305 @@
     first = last;
     --first;
-    size_t rank = n;
-    
+    typename traits::difference_type rank = n;
+
     double prev = 1.0;
     while (rank) {
@@ -286 +315 @@
       --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()));
   }
 
@@ -309 +355 @@
 
   template <class RandomAccessIterator>
-  double mad(RandomAccessIterator first, RandomAccessIterator last, 
+  double mad(RandomAccessIterator first, RandomAccessIterator last,
              bool sorted)
   {