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/concept_check.h"
00027 #include "yat/utility/DataWeight.h"
00028 #include "yat/utility/iterator_traits.h"
00029 #include "yat/utility/yat_assert.h"
00030 #include "yat/utility/WeightIterator.h"
00031
00032 #include <boost/concept_check.hpp>
00033
00034 #include <algorithm>
00035 #include <cmath>
00036 #include <limits>
00037 #include <numeric>
00038 #include <stdexcept>
00039 #include <vector>
00040
00041 namespace theplu {
00042 namespace yat {
00043 namespace statistics {
00044
00050 class Percentiler
00051 {
00052 public:
00060 Percentiler(double perc=50, bool sorted=false);
00061
00096 template<typename RandomAccessIterator>
00097 double operator()(RandomAccessIterator first,
00098 RandomAccessIterator last) const
00099 {
00100 BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
00101 BOOST_CONCEPT_ASSERT((utility::DataIteratorConcept<RandomAccessIterator>));
00102 if (first==last)
00103 return std::numeric_limits<double>::quiet_NaN();
00104 return calculate(first, last, sorted_,
00105 typename utility::weighted_iterator_traits<RandomAccessIterator>::type());
00106 }
00107 private:
00108 double perc_;
00109 bool sorted_;
00110
00111
00112 template<typename RandomAccessIterator>
00113 double calculate(RandomAccessIterator first, RandomAccessIterator last,
00114 bool sorted, utility::unweighted_iterator_tag tag) const;
00115
00116
00117 template<typename RandomAccessIterator>
00118 double calculate(RandomAccessIterator first, RandomAccessIterator last,
00119 bool sorted, utility::weighted_iterator_tag tag) const;
00120
00121
00122
00123
00124
00125 };
00126
00127
00128
00129
00130
00131 template<typename RandomAccessIterator>
00132 double
00133 Percentiler::calculate(RandomAccessIterator first,
00134 RandomAccessIterator last,
00135 bool sorted,
00136 utility::unweighted_iterator_tag tag) const
00137 {
00138
00139 if (first+1 == last)
00140 *first;
00141 if (sorted) {
00142 size_t n = last - first;
00143
00144 if (perc_>= 100.0)
00145 return *(--last);
00146 if (perc_<= 0.0)
00147 return *first;
00148 double j = n * perc_ / 100.0;
00149 size_t i = static_cast<size_t>(j);
00150 if (i==j)
00151 return (first[i]+first[i-1])/2;
00152 return first[i];
00153 }
00154
00155 std::vector<double> v_copy;
00156 v_copy.reserve(std::distance(first,last));
00157 std::copy(first, last, std::back_inserter(v_copy));
00158 std::sort(v_copy.begin(), v_copy.end());
00159 return calculate(v_copy.begin(), v_copy.end(), true, tag);
00160 }
00161
00162
00163
00164 template<typename RandomAccessIterator>
00165 double Percentiler::calculate(RandomAccessIterator first,
00166 RandomAccessIterator last,
00167 bool sorted,
00168 utility::weighted_iterator_tag tag) const
00169 {
00170 if (sorted) {
00171 utility::iterator_traits<RandomAccessIterator> trait;
00172 std::vector<double> accum_w;
00173 accum_w.reserve(last-first);
00174 std::partial_sum(weight_iterator(first),
00175 weight_iterator(last),
00176 std::back_inserter(accum_w));
00177
00178 double w_bound=perc_/100.0*accum_w.back();
00179 std::vector<double>::const_iterator upper(accum_w.begin());
00180 double margin=1e-10;
00181 while (upper!=accum_w.end() && *upper <= w_bound+margin)
00182 ++upper;
00183 while (upper!=accum_w.begin() &&
00184 (upper==accum_w.end() ||
00185 trait.weight(first+(upper-accum_w.begin()))==0.0))
00186 --upper;
00187 std::vector<double>::const_iterator lower(upper);
00188 while ( *(lower-1)>=w_bound-margin && lower>accum_w.begin())
00189 --lower;
00190
00191 return (trait.data(first+(upper-accum_w.begin()))+
00192 trait.data(first+(lower-accum_w.begin())))/2;
00193 }
00194
00195 std::vector<utility::DataWeight> v_copy;
00196 v_copy.reserve(last-first);
00197 std::copy(first, last, std::back_inserter(v_copy));
00198 std::sort(v_copy.begin(), v_copy.end());
00199 return calculate(v_copy.begin(), v_copy.end(), true, tag);
00200 }
00201
00202
00203 }}}
00204
00205 #endif