Changeset 1418 for trunk/yat/statistics


Ignore:
Timestamp:
Aug 20, 2008, 1:15:32 AM (13 years ago)
Author:
Peter
Message:

refs #366 - need to polish on some special cases

File:
1 edited

Legend:

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

    r1404 r1418  
    106106    // range is one value only is a special case
    107107    if (first+1 == last)
    108       return utility::iterator_traits<RandomAccessIterator>().data(first);
     108      *first;
    109109    if (sorted) {
    110110      size_t n = last - first;
    111111      // have to take care of this special case
    112       if (perc_>= 100.0 - 50.0/n)
     112      if (perc_>= 100.0)
    113113        return *(--last);
    114       if (perc_<= 50.0/n)
     114      if (perc_<= 0.0)
    115115        return *first;
    116       double j = n * (perc_-50.0/n) / 100.0;
    117       int i = static_cast<int>(j);
    118       return (1-j+floor(j))*first[i] + (j-floor(j))*first[i+1];
     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];
    119121    }
    120122
     
    142144                       std::back_inserter(accum_w));
    143145
    144       double sum_w=0;
    145146      double w_bound=perc_/100.0*accum_w.back();
    146       RandomAccessIterator upper(first);
    147       while (sum_w + 0.5*trait.weight(upper) < w_bound) {
    148         sum_w += trait.weight(upper);
     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)
    149150        ++upper;
    150       }
    151       if (upper==last)
    152         return trait.data(--last);
    153       if (sum_w==0.0) {
    154         while (first!=last)
    155           if (trait.weight(first))
    156             return trait.data(first);
    157         // all weights zero, hmm
    158         return 0;
    159       }
    160       // lower should have non-zero weight
    161       // and (lower, upper) should have zero weights
    162       RandomAccessIterator lower=upper-1;
    163       while (trait.weight(lower)==0)
     151      std::vector<double>::const_iterator lower(upper);
     152      while (lower!=accum_w.begin() && *lower > w_bound - allowed_error)
    164153        --lower;
    165 
    166       double p1=100.0/accum_w.back()*(sum_w + 0.5*trait.weight(upper));
    167       double p0=100.0/accum_w.back()*(sum_w - 0.5*trait.weight(lower));
    168       double alpha = (perc_-p0)/(p1-p0);
    169       utility::yat_assert<std::runtime_error>(alpha<=1.0,
    170                                               "Percentile: alpha out of bounds");
    171       return trait.data(lower)+alpha*(trait.data(upper)-trait.data(lower));
     154      return (trait.data(first+(upper-accum_w.begin()))+
     155              trait.data(first+(lower-accum_w.begin()+1)))/2;
    172156    }
    173157
Note: See TracChangeset for help on using the changeset viewer.