00001 #ifndef _theplu_yat_statistics_percentiler_
00002 #define _theplu_yat_statistics_percentiler_
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "yat/utility/DataWeight.h"
00027 #include "yat/utility/iterator_traits.h"
00028 #include "yat/utility/yat_assert.h"
00029 #include "yat/utility/WeightIterator.h"
00030
00031 #include <boost/concept_check.hpp>
00032
00033 #include <algorithm>
00034 #include <cmath>
00035 #include <numeric>
00036 #include <stdexcept>
00037 #include <vector>
00038
00039 namespace theplu {
00040 namespace yat {
00041 namespace statistics {
00042
00048 class Percentiler
00049 {
00050 public:
00058 Percentiler(double perc=50, bool sorted=false);
00059
00094 template<typename RandomAccessIterator>
00095 double operator()(RandomAccessIterator first,
00096 RandomAccessIterator last) const
00097 {
00098 BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
00099 return calculate(first, last, sorted_,
00100 typename utility::weighted_iterator_traits<RandomAccessIterator>::type());
00101 }
00102 private:
00103 double perc_;
00104 bool sorted_;
00105
00106
00107 template<typename RandomAccessIterator>
00108 double calculate(RandomAccessIterator first, RandomAccessIterator last,
00109 bool sorted, utility::unweighted_iterator_tag tag) const;
00110
00111
00112 template<typename RandomAccessIterator>
00113 double calculate(RandomAccessIterator first, RandomAccessIterator last,
00114 bool sorted, utility::weighted_iterator_tag tag) const;
00115
00116
00117
00118
00119
00120 };
00121
00122
00123
00124
00125
00126 template<typename RandomAccessIterator>
00127 double
00128 Percentiler::calculate(RandomAccessIterator first,
00129 RandomAccessIterator last,
00130 bool sorted,
00131 utility::unweighted_iterator_tag tag) const
00132 {
00133
00134 if (first+1 == last)
00135 *first;
00136 if (sorted) {
00137 size_t n = last - first;
00138
00139 if (perc_>= 100.0)
00140 return *(--last);
00141 if (perc_<= 0.0)
00142 return *first;
00143 double j = n * perc_ / 100.0;
00144 size_t i = static_cast<size_t>(j);
00145 if (i==j)
00146 return (first[i]+first[i-1])/2;
00147 return first[i];
00148 }
00149
00150 std::vector<double> v_copy;
00151 v_copy.reserve(std::distance(first,last));
00152 std::copy(first, last, std::back_inserter(v_copy));
00153 std::sort(v_copy.begin(), v_copy.end());
00154 return calculate(v_copy.begin(), v_copy.end(), true, tag);
00155 }
00156
00157
00158
00159 template<typename RandomAccessIterator>
00160 double Percentiler::calculate(RandomAccessIterator first,
00161 RandomAccessIterator last,
00162 bool sorted,
00163 utility::weighted_iterator_tag tag) const
00164 {
00165 if (sorted) {
00166 utility::iterator_traits<RandomAccessIterator> trait;
00167 std::vector<double> accum_w;
00168 accum_w.reserve(last-first);
00169 std::partial_sum(weight_iterator(first),
00170 weight_iterator(last),
00171 std::back_inserter(accum_w));
00172
00173 double w_bound=perc_/100.0*accum_w.back();
00174 std::vector<double>::const_iterator upper(accum_w.begin());
00175 double margin=1e-10;
00176 while (*upper <= w_bound+margin && upper!=accum_w.end())
00177 ++upper;
00178 if (upper==accum_w.end())
00179 --upper;
00180 std::vector<double>::const_iterator lower(upper);
00181 while ( *(lower-1)>=w_bound-margin && lower>accum_w.begin())
00182 --lower;
00183
00184 return (trait.data(first+(upper-accum_w.begin()))+
00185 trait.data(first+(lower-accum_w.begin())))/2;
00186 }
00187
00188 std::vector<utility::DataWeight> v_copy;
00189 v_copy.reserve(last-first);
00190 std::copy(first, last, std::back_inserter(v_copy));
00191 std::sort(v_copy.begin(), v_copy.end());
00192 return calculate(v_copy.begin(), v_copy.end(), true, tag);
00193 }
00194
00195
00196 }}}
00197
00198 #endif