Statistics/interface/lms_1d.h

0001 #ifndef CommonTools_Statistics_lms_1d_h
0002 #define CommonTools_Statistics_lms_1d_h
0003
0004 #include "CommonTools/Statistics/interface/StatisticsException.h"
0005
0006 #include <vector>
0007 #include <algorithm>
0008 #include <cmath>
0009
0010 /**
0011  *  Least Median Sum of Squares in one dimension.
0012  *  T must be sortable, 'addable', and dividable by 2.
0013  */
0014 template <class T>
0015 T lms_1d(std::vector<T> values) {
0016   switch (values.size()) {
0017     case 0:
0018       throw StatisticsException("Lms of empty vector undefined");
0019     case 1:
0020       return (T) * (values.begin());
0021     case 2:
0022       return (T)((*(values.begin()) + *(values.end() - 1)) / 2);
0023   };
0024   const int size = values.size();
0025   const int half = size / 2;
0026   const int n_steps = (size + 1) / 2;
0027
0028   sort(values.begin(), values.end());
0029
0030   T i_begin = *(values.begin());
0031   T i_end = *(values.begin() + half);
0032   T div = (T)fabs(i_begin - i_end);
0033   for (typename std::vector<T>::const_iterator i = values.begin(); i != (values.begin() + n_steps); ++i) {
0034     if (fabs((*i) - (*(i + half))) < div) {
0035       i_begin = (*i);
0036       i_end = *(i + half);
0037       div = (T)fabs(i_begin - i_end);
0038     }
0039   }
0040
0041   return (T)((i_begin + i_end) / 2);
0042 }
0043
0044 #endif