From e90a3b601bf90aa324e286735ac6944da6ef0dbe Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Fri, 10 Sep 2021 00:09:26 +0200
Subject: [PATCH] added the repeated median estimator of Siegel

---
 README.md          |  4 +++
 repeated_median.hh | 65 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 69 insertions(+)
 create mode 100644 repeated_median.hh

diff --git a/README.md b/README.md
index d0ca5ae..5fb1838 100644
--- a/README.md
+++ b/README.md
@@ -187,6 +187,10 @@ Implemented [Simple linear regression](https://en.wikipedia.org/wiki/Simple_line
 
 The [Theil-Sen estimator](https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator) is a [robust](https://en.wikipedia.org/wiki/Robust_statistics) [line fitting](https://en.wikipedia.org/wiki/Line_fitting) algorithm.
 
+### [repeated_median.hh](repeated_median.hh)
+
+The [repeated median estimator](https://en.wikipedia.org/wiki/Repeated_median_regression) is a [robust](https://en.wikipedia.org/wiki/Robust_statistics) [line fitting](https://en.wikipedia.org/wiki/Line_fitting) algorithm.
+
 ### [complex.hh](complex.hh)
 
 Faster alternative (no Inf/NaN handling) to the std::complex implementation.
diff --git a/repeated_median.hh b/repeated_median.hh
new file mode 100644
index 0000000..cc09cba
--- /dev/null
+++ b/repeated_median.hh
@@ -0,0 +1,65 @@
+/*
+Repeated median estimator of Siegel
+
+Copyright 2021 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include <algorithm>
+
+namespace DSP {
+
+template <typename TYPE, int SIZE>
+class RepeatedMedianEstimator
+{
+	TYPE inner_[SIZE-1], outer_[SIZE];
+	TYPE xint_, yint_, slope_;
+public:
+	RepeatedMedianEstimator() : 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) {
+			int count = 0;
+			for (int j = 0; j < LEN; ++j)
+				if (x[j] != x[i])
+					inner_[count++] = (x[j]*y[i] - x[i]*y[j]) / (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);
+		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;
+	}
+};
+
+}
+