Changeset 3236


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

BH correction on unsorted range. closes #796

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/statistics.cc

    r3137 r3236  
    4848using namespace theplu::yat;
    4949void test_benjamini_hochberg(test::Suite&);
     50void test_benjamini_hochberg_unsorted(test::Suite&);
    5051void test_entropy(test::Suite&);
    5152void test_mad(test::Suite&);
     
    135136  }
    136137  test_benjamini_hochberg(suite);
     138  test_benjamini_hochberg_unsorted(suite);
    137139  test_entropy(suite);
    138140  test_median_empty(suite);
     
    168170
    169171
     172void test_benjamini_hochberg_unsorted(test::Suite& suite)
     173{
     174  std::vector<double> p;
     175  p.push_back(0.015);
     176  p.push_back(0.0001);
     177  p.push_back(0.01);
     178  p.push_back(0.5);
     179  p.push_back(0.99);
     180  std::vector<double> q(p.size());
     181  statistics::benjamini_hochberg_unsorted(p.begin(), p.end(), q.begin());
     182  suite.add(suite.equal(q[1], p[1]*5));
     183  suite.add(suite.equal(q[2], p[2]*2.5));
     184  suite.add(suite.equal(q[0], 0.025));
     185  suite.add(suite.equal(q[3], p[3]*1.25));
     186  suite.add(suite.equal(q[4], 0.99));
     187
     188  // do nut run compiler test
     189  if (false) {
     190    using statistics::benjamini_hochberg_unsorted;
     191    boost::random_access_iterator_archetype<double> input;
     192    boost::mutable_random_access_iterator_archetype<double> result;
     193    benjamini_hochberg_unsorted(input, input, result);
     194  }
     195}
     196
     197
    170198void test_entropy(test::Suite& suite)
    171199{
  • 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.