From 3f3572f3774fd9185d9f46613c7d328698804ba4 Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Fri, 3 Sep 2021 01:22:02 +0200
Subject: [PATCH] added faster Theil-Sen estimator with random sampling

---
 theil_sen.hh | 46 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 46 insertions(+)

diff --git a/theil_sen.hh b/theil_sen.hh
index c565153..e6ae4bf 100644
--- a/theil_sen.hh
+++ b/theil_sen.hh
@@ -6,6 +6,8 @@ Copyright 2021 Ahmet Inan <inan@aicodix.de>
 
 #pragma once
 
+#include <random>
+#include <functional>
 #include <algorithm>
 
 namespace DSP {
@@ -52,5 +54,49 @@ public:
 	}
 };
 
+template <typename TYPE, int SIZE>
+class TheilSenEstimator2
+{
+	TYPE temp_[SIZE];
+	TYPE xint_, yint_, slope_;
+public:
+	TheilSenEstimator2() : xint_(0), yint_(0), slope_(0) {}
+	void compute(TYPE *x, TYPE *y, int LEN)
+	{
+		std::random_device rd;
+		std::default_random_engine generator(rd());
+		typedef std::uniform_int_distribution<int> uniform;
+		auto rand = std::bind(uniform(0, LEN-1), generator);
+		for (int i, j, k = 0; k < SIZE; ++k) {
+			while (x[i=rand()] == x[j=rand()]);
+			temp_[k] = (y[j] - y[i]) / (x[j] - x[i]);
+		}
+		std::nth_element(temp_, temp_+SIZE/2, temp_+SIZE);
+		slope_ = temp_[SIZE/2];
+		int count = std::min(LEN, SIZE);
+		for (int i = 0; i < count; ++i)
+			temp_[i] = y[i] - slope_ * x[i];
+		std::nth_element(temp_, temp_+count/2, temp_+count);
+		yint_ = temp_[count/2];
+		xint_ = - yint_ / slope_;
+	}
+	TYPE xint()
+	{
+		return xint_;
+	}
+	TYPE slope()
+	{
+		return slope_;
+	}
+	TYPE yint()
+	{
+		return yint_;
+	}
+	TYPE operator () (TYPE x)
+	{
+		return yint_ + slope_ * x;
+	}
+};
+
 }