Ignore:
Timestamp:
May 23, 2014, 3:42:51 PM (9 years ago)
Author:
Peter
Message:

BH correction on unsorted range. closes #796

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/statistics/utility.h

    r3137 r3236  
    99  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
    1010  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
    11   Copyright (C) 2009, 2010, 2011, 2013 Peter Johansson
     11  Copyright (C) 2009, 2010, 2011, 2013, 2014 Peter Johansson
    1212
    1313  This file is part of the yat library, http://dev.thep.lu.se/yat
     
    3636#include "yat/utility/deprecate.h"
    3737#include "yat/utility/iterator_traits.h"
     38#include "yat/utility/sort_index.h"
    3839#include "yat/utility/Vector.h"
    3940#include "yat/utility/VectorBase.h"
     
    4142
    4243#include <boost/concept_check.hpp>
     44#include <boost/iterator/permutation_iterator.hpp>
     45
    4346#include <gsl/gsl_math.h>
    4447#include <gsl/gsl_statistics_double.h>
     
    7275     constraint that \f$ q_m = p_m \f$.
    7376
    74      Requirements: \c BidirectionalIterator1 should be a
    75      \bidirectional_iterator and \c BidirectionalIterator2 should be a
    76      mutable \bidirectional_iterator
     77     Type Requirements:
     78     - \c BidirectionalIterator1 is a \bidirectional_iterator
     79     - \c BidirectionalIterator1::value type is convertible to \c double
     80     - \c BidirectionalIterator2 is a mutable \bidirectional_iterator
     81     - \c double is convertible to \c BidirectionalIterator2::value_type
    7782
    7883     \since New in yat 0.8
     
    8287                          BidirectionalIterator1 last,
    8388                          BidirectionalIterator2 result);
     89
     90
     91  /**
     92     \brief Benjamini Hochberg multiple test correction
     93
     94     Similar to benjamini_hochberg() but does not assume that input
     95     range, [first, last), is sorted. The resulting range is the same
     96     as if sorting input range, call benjamini_hochberg, and unsort
     97     the result range.
     98
     99     Type Requirements:
     100     - \c RandomAccessIterator models a \random_access_iterator
     101     - \c RandomAccessIterator::value type is convertible to \c double
     102     - \c MutableRandomAccessIterator models a mutable \random_access_iterator
     103     - \c double is convertible to \c MutableRandomAccessIterator::value_type
     104
     105     \since New in yat 0.13
     106   */
     107  template<typename RandomAccessIterator, typename MutableRandomAccessIterator>
     108  void benjamini_hochberg_unsorted(RandomAccessIterator first,
     109                                   RandomAccessIterator last,
     110                                   MutableRandomAccessIterator result);
    84111
    85112
     
    95122  /// @return cumulative hypergeomtric distribution functions P(k).
    96123  ///
    97   double cdf_hypergeometric_P(unsigned int k, unsigned int n1, 
     124  double cdf_hypergeometric_P(unsigned int k, unsigned int n1,
    98125                              unsigned int n2, unsigned int t);
    99126
     
    248275  ///
    249276  double skewness(const utility::VectorBase&);
    250  
     277
    251278
    252279  //////// template implementations ///////////
     
    259286    for (size_t i=0; first!=last; ++i, ++first)
    260287      o.add(traits.data(first), target.binary(i), traits.weight(first));
    261   } 
     288  }
    262289
    263290
     
    268295  {
    269296    using boost::Mutable_BidirectionalIterator;
    270     BOOST_CONCEPT_ASSERT((boost::BidirectionalIterator<BidirectionalIterator1>));
     297    using boost::BidirectionalIterator;
     298    BOOST_CONCEPT_ASSERT((BidirectionalIterator<BidirectionalIterator1>));
    271299    BOOST_CONCEPT_ASSERT((Mutable_BidirectionalIterator<BidirectionalIterator2>));
    272     size_t n = std::distance(first, last);
     300    typedef typename std::iterator_traits<BidirectionalIterator1> traits;
     301    typename traits::difference_type n = std::distance(first, last);
    273302    if (!n)
    274303      return;
     
    276305    first = last;
    277306    --first;
    278     size_t rank = n;
    279    
     307    typename traits::difference_type rank = n;
     308
    280309    double prev = 1.0;
    281310    while (rank) {
     
    286315      --result;
    287316    }
     317  }
     318
     319
     320  template<typename RandomAccessIterator, typename MutableRandomAccessIterator>
     321  void benjamini_hochberg_unsorted(RandomAccessIterator first,
     322                                   RandomAccessIterator last,
     323                                   MutableRandomAccessIterator result)
     324  {
     325    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
     326    using boost::Mutable_RandomAccessIterator;
     327    BOOST_CONCEPT_ASSERT((Mutable_RandomAccessIterator<MutableRandomAccessIterator>));
     328
     329    std::vector<size_t> idx;
     330    utility::sort_index(first, last, idx);
     331    benjamini_hochberg(boost::make_permutation_iterator(first, idx.begin()),
     332                       boost::make_permutation_iterator(first, idx.end()),
     333                       boost::make_permutation_iterator(result, idx.begin()));
    288334  }
    289335
     
    309355
    310356  template <class RandomAccessIterator>
    311   double mad(RandomAccessIterator first, RandomAccessIterator last, 
     357  double mad(RandomAccessIterator first, RandomAccessIterator last,
    312358             bool sorted)
    313359  {
Note: See TracChangeset for help on using the changeset viewer.