source: trunk/yat/statistics/Percentiler.h @ 1418

Last change on this file since 1418 was 1418, checked in by Peter, 13 years ago

refs #366 - need to polish on some special cases

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 4.8 KB
Line 
1#ifndef _theplu_yat_statistics_percentiler_
2#define _theplu_yat_statistics_percentiler_
3
4// $Id: Percentiler.h 1418 2008-08-19 23:15:32Z peter $
5
6/*
7  Copyright (C) 2008 Peter Johansson
8
9  This file is part of the yat library, http://trac.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "yat/utility/DataWeight.h"
28#include "yat/utility/iterator_traits.h"
29#include "yat/utility/yat_assert.h"
30#include "yat/utility/WeightIterator.h"
31
32#include <algorithm>
33#include <cassert>
34#include <cmath>
35#include <numeric>
36#include <stdexcept>
37#include <vector>
38
39#include <iostream>
40
41namespace theplu {
42namespace yat {
43namespace statistics { 
44
45  /**
46     \brief Functor to calculate percentile of a range
47
48     \since New in yat 0.5
49   */
50  class Percentiler
51  {
52  public:
53    /**
54       \param perc percentile to calculate [0,100]. Defaule value is
55       50, which implies class will calculate median.
56       \param sorted if true class assumes that ranges are already
57       sorted, if false the range will copied to a new range which is
58       sorted.
59     */
60    Percentiler(double perc=50, bool sorted=false);
61
62    /**
63       Function is a non-mutable function, i.e., \a first and \a last
64       can be const_iterators.
65
66       \return percentile of range
67     */
68    template<typename RandomAccessIterator>
69    double operator()(RandomAccessIterator first, 
70                      RandomAccessIterator last) const
71    {
72      return calculate(first, last, sorted_, 
73       typename utility::weighted_iterator_traits<RandomAccessIterator>::type());
74    }
75  private:
76    double perc_;
77    bool sorted_;
78
79    // unweighted version
80    template<typename RandomAccessIterator>
81    double calculate(RandomAccessIterator first, RandomAccessIterator last,
82                     bool sorted, utility::unweighted_iterator_tag tag) const;
83
84    // weighted version
85    template<typename RandomAccessIterator>
86    double calculate(RandomAccessIterator first, RandomAccessIterator last,
87                     bool sorted, utility::weighted_iterator_tag tag) const;
88
89    // using compiler generated copy
90    //Percentiler(const Percentiler&);
91    //Percentiler& operator=(const Percentiler&);
92
93  };
94
95  // template implementations
96  // /////////////////////////
97
98  // unweighted version
99  template<typename RandomAccessIterator>
100  double 
101  Percentiler::calculate(RandomAccessIterator first, 
102                         RandomAccessIterator last,
103                         bool sorted,
104                         utility::unweighted_iterator_tag tag) const
105  {
106    // range is one value only is a special case
107    if (first+1 == last)
108      *first;
109    if (sorted) {
110      size_t n = last - first;
111      // have to take care of this special case
112      if (perc_>= 100.0)
113        return *(--last);
114      if (perc_<= 0.0)
115        return *first;
116      double j = n * perc_ / 100.0;
117      size_t i = static_cast<size_t>(j);
118      if (i==j)
119        return (first[i]+first[i-1])/2;
120      return first[i];
121    }
122
123    std::vector<double> v_copy;
124    v_copy.reserve(std::distance(first,last));
125    std::copy(first, last, std::back_inserter(v_copy));
126    std::sort(v_copy.begin(), v_copy.end());
127    return calculate(v_copy.begin(), v_copy.end(), true, tag);
128  }
129
130
131  // weighted version
132  template<typename RandomAccessIterator>
133  double Percentiler::calculate(RandomAccessIterator first, 
134                                RandomAccessIterator last,
135                                bool sorted,
136                                utility::weighted_iterator_tag tag) const
137  {
138    if (sorted) {
139      utility::iterator_traits<RandomAccessIterator> trait;
140      std::vector<double> accum_w;
141      accum_w.reserve(last-first);
142      std::partial_sum(weight_iterator(first),
143                       weight_iterator(last),
144                       std::back_inserter(accum_w));
145
146      double w_bound=perc_/100.0*accum_w.back();
147      std::vector<double>::const_iterator upper(accum_w.begin());
148      double allowed_error=1e-10;
149      while (upper != accum_w.end() && *upper < w_bound + allowed_error)
150        ++upper;
151      std::vector<double>::const_iterator lower(upper);
152      while (lower!=accum_w.begin() && *lower > w_bound - allowed_error)
153        --lower;
154      return (trait.data(first+(upper-accum_w.begin()))+
155              trait.data(first+(lower-accum_w.begin()+1)))/2;
156    }
157
158    std::vector<utility::DataWeight> v_copy;
159    v_copy.reserve(last-first);
160    std::copy(first, last, std::back_inserter(v_copy));
161    std::sort(v_copy.begin(), v_copy.end());
162    return calculate(v_copy.begin(), v_copy.end(), true, tag);
163  }
164
165
166}}} // of namespace statistics, yat, and theplu
167
168#endif
Note: See TracBrowser for help on using the repository browser.