diff --git a/repeated_median.hh b/repeated_median.hh index cc09cba..2ba25b7 100644 --- a/repeated_median.hh +++ b/repeated_median.hh @@ -61,5 +61,50 @@ public: } }; +template +class RepeatedMedianEstimator2 +{ + TYPE inner_[SIZE-1], outer_[SIZE]; + TYPE xint_, yint_, slope_; +public: + RepeatedMedianEstimator2() : xint_(0), yint_(0), slope_(0) {} + void compute(TYPE *x, TYPE *y, int LEN) + { + if (LEN > SIZE) + LEN = SIZE; + for (int i = 0; i < LEN; ++i) { + int count = 0; + for (int j = 0; j < LEN; ++j) + if (x[j] != x[i]) + inner_[count++] = (y[j] - y[i]) / (x[j] - x[i]); + std::nth_element(inner_, inner_+count/2, inner_+count); + outer_[i] = inner_[count/2]; + } + std::nth_element(outer_, outer_+LEN/2, outer_+LEN); + slope_ = outer_[LEN/2]; + for (int i = 0; i < LEN; ++i) + outer_[i] = y[i] - slope_ * x[i]; + std::nth_element(outer_, outer_+LEN/2, outer_+LEN); + yint_ = outer_[LEN/2]; + xint_ = - yint_ / slope_; + } + TYPE xint() + { + return xint_; + } + TYPE slope() + { + return slope_; + } + TYPE yint() + { + return yint_; + } + TYPE operator () (TYPE x) + { + return yint_ + slope_ * x; + } +}; + }