通过泰尔森估计器在熊猫时间序列中绘制回归线的置信区间

Question

我正在开发一个时间序列数据框，例如：

df = pd.DataFrame({'year':[ '1990','1991','1992','1993','1994','1995','1996','1997','1998','1999','2000'],
                           'count':[96,184,148,154,160,149,124,274,322,301,300]})

我有兴趣使用Theil-Sen's slope estimator 和Scipy 中的函数找到回归线

https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.stats.mstats.theilslopes.html

如下codesn-p:

from scipy import stats                                         


x = df['year'].astype(float)
y = df['count']


res = stats.theilslopes(y, x, 0.90)


fig = plt.figure()
ax = fig.add_subplot(111)

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

ax.plot(x, res[1] + res[0] * x, 'r-')


plt.ylabel('count')
plt.xlabel('year')

plt.show()

我可以找到 Theil-Sen's regression line，但是在添加可视化 confidence interval 的代码时，可视化有些不现实。

    ax.plot(x, res[1] + res[2] * x, 'r--')
    ax.plot(x, res[1] + res[3] * x, 'r--')

有没有办法，我可以用theil-sen's slope estimator和confidence interval得到一个可视化，如下图：

我们将非常感谢您的帮助。 谢谢！

Answer 1

查看documentation，您必须获得自己的置信区间估计，因为没有给出任何置信区间估计，并且您不能仅基于斜率不确定性。误差传播或蒙特卡洛是您可以用来估计置信区间的一些方法。

但是，您的图看起来不正常的原因是归一化，如果您在拟合之前对 x 轴进行归一化，您会得到一个看起来正常的置信区间，而不是一个会爆炸的置信区间。然而，这不是一个真正的置信区间，因为唯一的不确定性是斜率，因此它在截距处交叉为零。在这里，我选择标准化为 -1,1 以生成您想要的图形。但我也提供了一些辅助函数来帮助您在未来以不同的方式执行此操作。

import pandas as pd
df = pd.DataFrame({'year':[ '1990','1991','1992','1993','1994','1995','1996','1997','1998','1999','2000'],
                           'count':[96,184,148,154,160,149,124,274,322,301,300]})

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

x = df['year'].astype(float)
y = df['count']

def norm0_1(x):
    return (x - x.min())/ (x.max()-x.min())

def denorm(x_norm, x):
    return (x_norm * (x.max() - x.min()) + x.min())

x_norm = 2*norm0_1(x)-1

res = stats.theilslopes(y, x_norm, 0.90)


fig = plt.figure()
ax = fig.add_subplot(111)

ax.plot(x, y, 'b.')
ax.plot(x, res[1] + res[0] * x_norm, 'r-')
ax.plot(x, res[1] + res[2] * x_norm, 'r--')
ax.plot(x, res[1] + res[3] * x_norm, 'r--')


plt.ylabel('count')
plt.xlabel('year')

plt.show()

标准版：

您的代码生成的默认版本：