Python SciPy stats.theilslopes用法及代碼示例

本文簡要介紹 python 語言中 scipy.stats.theilslopes 的用法。

用法:

scipy.stats.theilslopes(y, x=None, alpha=0.95, method='separate')#

計算一組點 (x, y) 的 Theil-Sen 估計量。

theilslopes 實現了穩健線性回歸的方法。它將斜率計算為配對值之間所有斜率的中值。

參數 ：：

y： array_like

因變量。

x： 數組 或 None，可選

自變量。如果沒有，請改用arange(len(y))。

alpha： 浮點數，可選

置信度介於 0 和 1 之間。默認為 95% 置信度。請注意， alpha 在 0.5 附近是對稱的，即 0.1 和 0.9 都被解釋為“找到 90% 置信區間”。

method： {‘joint’, ‘separate’}，可選

用於計算截距估計值的方法。支持以下方法，

• ‘joint’: Uses np.median(y - slope * x) as intercept.

• ‘separate’: Uses np.median(y) - slope * np.median(x)

  as intercept.

默認值為‘separate’。

返回 ：：

result： TheilslopesResult 實例

返回值是一個具有以下屬性的對象：

坡 浮點數

泰爾坡。

截距 浮點數

Theil線的攔截。

low_slope 浮點數

斜率置信區間的下限。

high_slope 浮點數

斜率置信區間的上限。

注意：

theilslopes 的實現遵循 [1]。截距在[1]中沒有定義，這裏定義為median(y) - slope*median(x)，在[3]中給出。截距的其他定義存在於文獻中，例如[4]中的median(y - slope*x)。計算截距的方法可以由參數 method 確定。沒有給出截距的置信區間，因為[1]中沒有解決這個問題。

為了與舊版本的 SciPy 兼容，返回值的行為類似於長度為 4 的 namedtuple ，其中包含字段 slope 、 intercept 、 low_slope 和 high_slope ，因此可以繼續編寫：

slope, intercept, low_slope, high_slope = theilslopes(y, x)

參考：

[1] (1,2,3)

PK Sen，“基於 Kendall tau 的回歸係數估計”，J. Am。統計。協會，卷。 63，第 1379-1389 頁，1968 年。

[2]

H. Theil，“線性和多項式回歸分析 I、II 和 III 的rank-invariant 方法”，Nederl。阿卡德。韋滕施，PROC。 53：第 386-392 頁、第 521-525 頁、第 1397-1412、1950 頁。

[3]

W.L. Conover，“Practical nonparametric statistics”，第 2 版，John Wiley and Sons，紐約，第 493 頁。

[4]

https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator

例子：

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

>>> x = np.linspace(-5, 5, num=150)
>>> y = x + np.random.normal(size=x.size)
>>> y[11:15] += 10  # add outliers
>>> y[-5:] -= 7

計算斜率、截距和 90% 置信區間。為了進行比較，還可以使用 linregress 計算最小二乘擬合：

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

繪製結果。 Theil-Sen 回歸線以紅色顯示，紅色虛線表示斜率的置信區間(請注意，紅色虛線不是回歸的置信區間，因為不包括截距的置信區間)。綠線顯示用於比較的最小二乘擬合。

>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> ax.plot(x, y, 'b.')
>>> ax.plot(x, res[1] + res[0] * x, 'r-')
>>> ax.plot(x, res[1] + res[2] * x, 'r--')
>>> ax.plot(x, res[1] + res[3] * x, 'r--')
>>> ax.plot(x, lsq_res[1] + lsq_res[0] * x, 'g-')
>>> plt.show()