时间序列预测：预测间隔

2023-02-04 13:26 474 阅读 ID：756
磐创AI
磐创AI
主页

在现实世界中进行预测是一项重要的任务。考虑预测能源需求、温度、食物供应和健康指标，仅举几个例子。在这些情况下获得不准确的预测会对人们的生活产生重大影响。

这就是预测间隔可以帮助的地方。预测区间用于提供预测可能具有特定置信度的范围。例如，如果你以95%的置信度进行了100个预测，那么100个预测中有95个在预测区间内。通过使用预测区间，你可以解释预测中的不确定性以及数据的随机变化。

加载温度数据

在本文中，我们将预测巴西圣保罗的月平均气温。数据集由NCEI收集和管理. 感谢所有参与的人与世界分享这个数据集！

不幸的是，在现实世界中，数据从来都不是你想要的格式。

在下面的单元格中，我们使用Pandas将其转换为两列。日期和温度。结果数据使用Matplotlib可视化。

# ---- Packages ----
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import normal_ad
from scipy import stats
from sklearn import metrics

!pip install pmdarima --quiet
import pmdarima as pm
pd.plotting.register_matplotlib_converters()

# ---- Data Transformations ----
months_dict = {'JAN':1,'FEB':2,'MAR':3,'APR':4,'MAY':5,'JUN':6,
               'JUL':7,'AUG':8,'SEP':9,'OCT':10,'NOV':11,'DEC':12}

df = pd.read_csv('../input/temperature-timeseries-for-some-brazilian-cities/station_sao_paulo.csv')
df = df[['YEAR'] + list(months_dict.keys())]

df_sp = pd.melt(df, 
        id_vars='YEAR',
        value_vars=months_dict.keys(),
        var_name='MONTH', 
        value_name='Sum of Value').replace({"MONTH": months_dict})

df_sp['DAY'] = 1 #need a day for the pd.to_datetime() function
df_sp['DATE'] = pd.to_datetime(df_sp[['YEAR','MONTH','DAY']])
df_sp = df_sp[['DATE','Sum of Value']].rename(columns={'DATE':'date','Sum of Value':'temp'})
df_sp = df_sp.sort_values('date').reset_index(drop=True)

# ---- Visualize Data ----
plt.figure(figsize=(15,7))
plt.title("Sao Paulo AVG Monthly Temperature - w/ Err")
plt.plot(df_sp["date"], df_sp['temp'], color='#379BDB', label='Original')
plt.show(block=False)

在上面的图中，我们可以看到数据集缺少值和错误。

我们可以简单地对~1980到~2015之间的372个数据点进行采样，这些数据点没有误差。我们将在本文的其余部分使用此数据示例。

df_sp = df_sp[(df_sp['date'] > '1982-03-01') & (df_sp['date'] < '2013-04-01')]
df_sp = df_sp.set_index(df_sp['date'], drop=True).drop(columns=["date"])

train, valid = df_sp[:int(len(df_sp)*0.8)], df_sp[int(len(df_sp)*0.8):]

plt.figure(figsize=(15,7))
plt.title("Sao Paulo AVG Monthly Temperature - 1982-03-01 to 2013-04-01")
plt.plot(train.index, train['temp'], color='#379BDB', label='Train')
plt.plot(valid.index, valid['temp'], color='#fc7d0b', label='Valid')
plt.xlabel('Date')
plt.ylabel('Temp - Celsius')
plt.legend()
plt.show(block=False)

创建模型

我使用过许多时间序列模型，但我总是回到ARIMA模型。这些模型是可靠的，并且经常优于竞争模型类型（NeuralProphet、ExponentialSmoothing、Last Value）。

我们将使用ARIMA的修改版本，称为SARIMA。SARIMA在ARIMA中添加了一个滞后项，用于跟踪数据的季节性。

使用名为pmdarima的包⁹ 我们可以自动调整模型参数。有关ARIMA模型的更详细解释，请查看本文。

https://towardsdatascience.com/time-series-forecasting-with-arima-sarima-and-sarimax-ee61099e78f6

SARIMA_model = pm.auto_arima(train['temp'], start_p=1, start_q=1,
                         test='adf',
                         max_p=3, max_q=3, 
                         m=12, #annual frequency(12 for month, 7 for week etc) 
                         start_P=0, 
                         seasonal=True, #set to seasonal
                         d=None, 
                         D=1, #order of the seasonal differencing
                         trace=False,
                         error_action='ignore',  
                         suppress_warnings=True, 
                         stepwise=True)

上面的单元格为我们提供了适合ARIMA模型的最佳顺序和季节顺序。

在下面的单元中，我们只做了这一点，并对验证数据集迭代地进行预测。

sarima_preds = []
df_sp.index = df_sp.index.to_period('M') 

for i in range(len(valid)):
    m_sarima = ARIMA(df_sp[:len(train)+i]['temp'], order=(1,0,0), seasonal_order=(0, 1, 2, 12)).fit()
    sarima_preds.append(m_sarima.forecast(1).values[0])

residuals = sorted([x - y for x, y in zip(sarima_preds, valid['temp'].values)])

预测间隔

方法1:RMSFE

我们可以使用的第一种方法叫做RMSFE（均方根预测误差）。RMSFE与RMSE非常相似。唯一的区别是，RMSFE必须根据对未知数据（即验证或测试集）的预测的残差项进行计算。

需要注意的是，只有假设验证预测的残差是正态分布的，我们才能使用这种方法。为了确定情况是否如此，我们将使用PP图，并用K方检验检验其正态性。

PP图根据正态分布图绘制数据样本，如果正态分布，数据点将形成一条直线。

三个正态性检验使用p值确定数据样本来自正态分布总体的可能性。每个测试的无效假设是“样本来自正态分布”。这意味着，如果得到的p值低于选定的α值，则否定空假设。因此，有证据表明数据来自非正态分布。对于本文，我们将使用0.01的Alpha值。

sw_result = stats.shapiro(residuals)
ad_result = normal_ad(np.array(residuals), axis=0)
dag_result = stats.normaltest(residuals, axis=0, nan_policy='propagate')

plt.figure(figsize=(15,7))
res = stats.probplot(residuals, plot=plt)
ax = plt.gca()
ax.annotate("SW p-val: {:.4f}".format(sw_result[1]), xy=(0.05,0.9), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))
ax.annotate("AD p-val: {:.4f}".format(ad_result[1]), xy=(0.05,0.8), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))
ax.annotate("DAG p-val: {:.4f}".format(dag_result[1]), xy=(0.05,0.7), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))

plt.show()

太棒了所有三个测试都返回了大于0.01的α值的p值。这意味着不能拒绝零假设，数据点很可能来自正态分布。我们现在可以使用RMSFE生成预测区间。

这里的第一步是选择我们想要提供的信心程度。我们希望我们的预测在75%、95%或99%的预测区间内吗？我们将使用95%的预测区间。

在正态分布中，95%的数据点落在平均值的1.96个标准差内，因此我们将1.96乘以RMSFE以获得预测区间大小。如下图所示。

RMSFE = np.sqrt(sum([x**2 for x in residuals]) / len(residuals))
band_size = 1.96*RMSFE

fig, ax = plt.subplots(figsize=(15,7))
ax.plot(valid.index, valid['temp'], color='#fc7d0b', label='Valid')
ax.scatter(valid.index, sarima_preds)
ax.fill_between(valid.index, (valid['temp']-band_size), (valid['temp']+band_size), color='b', alpha=.1)
ax.set_title("Predictions w/ 95% Confidence")
ax.set_xlabel('Date')
ax.set_ylabel('Temp - Celsius')
plt.show()

该方法的缺点是预测区间高度依赖于验证预测的残差。

方法2:BCVR

注：BCVR是我提出的一种理论方法！

让我们将此方法称为BCVR（Bootstrapping Cross Validation Residuals）。BCVR试图从交叉验证和Bootstrap中获益。Bootstrapping残差是生成预测区间的常用方法，通常产生与具有正态分布残差的RMSFE方法相似的结果，但在非正态残差上的表现略优于RMSFE。此外，通过使用交叉验证，BCVR应生成更能代表整个数据集的残差分布。

我们可以通过执行交叉验证来生成残差。我们随机选择一个长度在250到372点之间的训练样本，并进行一步预测。然后，我们计算该预测的残差，并重复该过程1000次。

crossval_count = 1000
min_size = 250
max_random = len(df_sp) - (min_size)

crossval_resids = []

for b in range(crossval_count):
    start, end = [min_size*i+x for i, x in enumerate(sorted([np.random.randint(max_random) for y in range(2)]))]

    m_sarima = ARIMA(df_sp[start:end]['temp'], order=(1,0,0), seasonal_order=(0, 1, 2, 12), enforce_invertibility=False).fit()
    pred = m_sarima.forecast(1).values[0]
    crossval_resids.append(pred - df_sp[end:end+1]['temp'].values[0])

bsed_residuals = sorted(crossval_resids)

然后，我们可以测试残差分布的正态性。

sw_result = stats.shapiro(bsed_residuals)
ad_result = normal_ad(np.array(bsed_residuals), axis=0)
dag_result = stats.normaltest(bsed_residuals, axis=0, nan_policy='propagate')

plt.figure(figsize=(15,7))
res = stats.probplot(bsed_residuals, plot=plt)
ax = plt.gca()
ax.annotate("SW p-val: {:.4f}".format(sw_result[1]), xy=(0.05,0.9), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))
ax.annotate("AD p-val: {:.4f}".format(ad_result[1]), xy=(0.05,0.8), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))
ax.annotate("DAG p-val: {:.4f}".format(dag_result[1]), xy=(0.05,0.7), xycoords='axes fraction', fontsize=15,
            bbox=dict(boxstyle="round", fc="none", ec="gray", pad=0.6))

plt.show()

哦，不！我们的p值远远低于Alpha阈值，所以应该拒绝零假设吗？不用担心，我们可以使用一种称为bootstrap的技术来测量剩余方差。为此，我们采用交叉验证的残差，并使用替换进行随机抽样。然后，我们计算重采样集的标准差，并将其存储在数组中。我们重复这个过程几次，然后取存储的bootstrap标准差的平均值/中值。

我在这个实现中使用提前停止来减少计算需求。在下面的单元中，我们每200次迭代检查一次引导样本的平均标准差。如果平均标准偏差与先前测量值的偏差不超过0.001%，则我们终止循环。

在我们的例子中，停止发生在第7500次迭代上。

bs_stds=[]
last_std = None

for i in range(10000):
    samples = [bsed_residuals[np.random.randint(len(bsed_residuals))] for x in range(len(bsed_residuals))]
    bs_stds.append(np.std(samples))
    if i%500 == 0 and i!=0:
        if i == 500:
            last_std = np.median(bs_stds)
        else:
            current_std = np.median(bs_stds)
            if abs(current_std - last_std) < last_std*0.00001:
                print("Early Stopping Iteration: {}".format(i))
                break
            else:
                last_std = np.median(bs_stds)        

bs_std = np.median(bs_stds)
print("Median of Bootstrapped STD's: {:.4f}".format(bs_std))

bs_band_size = 1.96*bs_std 

fig, ax = plt.subplots(figsize=(15,7))
ax.plot(valid.index, valid['temp'], color='#fc7d0b', label='Valid')
ax.scatter(valid.index, sarima_preds)
ax.fill_between(valid.index, (valid['temp']-bs_band_size), (valid['temp']+bs_band_size), color='b', alpha=.1)
ax.set_title("Predictions w/ 95% Confidence")
ax.set_xlabel('Date')
ax.set_ylabel('Temp - Celsius')
plt.show()

BCVR方法的缺点是一台机器的计算成本很高。谢天谢地，这种方法非常适合集群计算。这意味着可以在多台机器之间分配工作负载。

例如，如果我们有一个有10个节点的集群，并且想要执行1000个引导样本，我们可以让每个节点同时执行100个样本。这将大大减少计算时间，并允许我们增加bootstrap样本的数量。

RMSFE与BCVR

这两种方法比较如何？

bs_band_size = 1.96*bs_std 

fig, ax = plt.subplots(figsize=(15,7))
ax.fill_between(valid.index, (valid['temp']-band_size), (valid['temp']+band_size), color='blue', alpha=0.1)
ax.fill_between(valid.index, (valid['temp']-bs_band_size), (valid['temp']+bs_band_size), color='yellow', alpha=0.3)

ax.plot(valid.index, valid['temp'], color='#fc7d0b', label='Valid')
ax.scatter(valid.index, sarima_preds)

ax.set_title("Predictions w/ 95% Confidence")
ax.set_xlabel('Date')
ax.set_ylabel('Temp - Celsius')

patches = [Patch(color='blue', alpha=0.1), Patch(color='#ececb2')]
labels = ['RMSFE','Bootstrapping Cross-Validation Residuals']
ax.legend(patches, labels)

plt.show()

在这两种情况下，所有75个点都在预测区间内。虽然不太可能（每种可能性为2%），但仍有一些点非常接近边界。这两种方法的间隔范围略有不同。

BCVR方法产生的预测区间范围比RMSFE方法略窄（参见上面围绕BCVR预测区间的淡蓝色线）。

最后

该示例表明，所提出的BCVR实现可能有助于生成更能代表整个样本空间的预测区间。看到它在不同场景中的应用会很有趣。

本文的代码可以在这里找到。

https://github.com/brendanartley/Medium-Article-Code/blob/main/code/time-series-forecasting-prediction-intervals.ipynb

参考引用

  • Forecasting: Principles and Practice — Rob J Hyndman and George Athanasopoulos
  • Confidence and Prediction intervals for forecasted values — Charles Zaiontz
  • Add Prediction Intervals to your Forecasting Model — Marco Cerliani
  • Bootstrapping Prediction Intervals — Dan Saattrup Nielson
  • A Gentle Intro to Normality Tests — Jason Brownlee
  • NCEI (National Centers for Environmental Information)
  • NASA GISS (Goddard Institute for Space Studies)
  • Shapiro-Wilk Test — Wikipedia
  • Pmdarina

感谢阅读！

 收起内容
免责声明：作者保留权利，不代表本站立场。如想了解更多和作者有关的信息可以查看页面右侧作者信息卡片。
    快抢沙发！
      加载更多评论
    反馈
    to-top--btn