https://stats.stackexchange.com/a/366/69508

366_69508.py

from scipy import stats
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(-1, 1, num=150)

y = 10*x + np.random.normal(size=x.size)

y[11:15] += 10
y[-5:] -= 7

res = stats.theilslopes(y, x, 0.90, method='separate')
lsq_res = stats.linregress(x, y)
fig = plt.figure()

ax = fig.add_subplot(111)

ax.plot(x, y, 'k.')

ax.plot(x, res[1] + res[0] * x, 'y-', label='Theil-Sen')

ax.plot(x, lsq_res[1] + lsq_res[0] * x, 'm-', label='Least Squares')

plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()