• 【ML】使用 ARIMA 进行时间序列预测

     🔎大家好，我是Sonhhxg_柒，希望你看完之后，能对你有所帮助，不足请指正！共同学习交流🔎

    📝个人主页－Sonhhxg_柒的博客_CSDN博客 📃

    🎁欢迎各位→点赞👍 + 收藏⭐️ + 留言📝​

    📣系列专栏 - 机器学习【ML】 自然语言处理【NLP】  深度学习【DL】

    ​

     🖍foreword

    ✔说明⇢本人讲解主要包括Python、机器学习（ML）、深度学习（DL）、自然语言处理（NLP）等内容。

    如果你对这个系列感兴趣的话，可以关注订阅哟👋

    介绍

    在本课中，您将发现一种使用ARIMA构建模型的特定方法：自动回归集成移动平均。ARIMA 模型特别适合拟合显示非平稳性的数据。

    一般概念

    为了能够使用 ARIMA，您需要了解一些概念：

    • 🎓 平稳性。从统计学的角度来看，平稳性是指数据的分布在随时间变化时不会发生变化。因此，非平稳数据会显示由于趋势而导致的波动，必须对其进行转换才能进行分析。例如，季节性可以引入数据波动，并且可以通过“季节性差异”过程消除。

    • 🎓 差异化。再次从统计上下文来看，差分数据是指通过消除其非恒定趋势来转换非平稳数据以使其平稳的过程。“差分消除了时间序列水平的变化，消除了趋势和季节性，从而稳定了时间序列的平均值。” 士雄等人的论文

    时间序列背景下的 ARIMA

    让我们拆开 ARIMA 的各个部分，以更好地了解它如何帮助我们对时间序列进行建模并帮助我们对其进行预测。

    • AR - 自回归。顾名思义，自回归模型会及时“回溯”以分析数据中的先前值并对它们做出假设。这些先前的值称为“滞后”。一个例子是显示铅笔每月销售额的数据。每个月的销售总额将被视为数据集中的“不断变化的变量”。这个模型被构建为“感兴趣的不断变化的变量根据其自身的滞后（即先验）值进行回归”。维基百科

    • I-用于集成。与类似的“ARMA”模型相反，ARIMA 中的“I”指的是其集成方面。当应用差分步骤以消除非平稳性时，数据是“集成的”。

    • MA - 移动平均线。该模型的移动平均方面是指通过观察当前和过去的滞后值确定的输出变量。

    底线：ARIMA 用于使模型尽可能接近时间序列数据的特殊形式。

    练习 - 构建 ARIMA 模型

    打开本课中的/working文件夹，找到notebook.ipynb文件。

    1. 运行 notebook 以加载statsmodelsPython 库；对于 ARIMA 模型，您将需要它。

    2. 加载必要的库

    3. 现在，加载更多用于绘制数据的库：

      1. import os
      2. import warnings
      3. import matplotlib.pyplot as plt
      4. import numpy as np
      5. import pandas as pd
      6. import datetime as dt
      7. import math
      8.  
      9. from pandas.plotting import autocorrelation_plot
      10. from statsmodels.tsa.statespace.sarimax import SARIMAX
      11. from sklearn.preprocessing import MinMaxScaler
      12. from common.utils import load_data, mape
      13. from IPython.display import Image
      14.  
      15. %matplotlib inline
      16. pd.options.display.float_format = '{:,.2f}'.format
      17. np.set_printoptions(precision=2)
      18. warnings.filterwarnings("ignore") # specify to ignore warning messages

    4. 将文件中的数据加载/data/energy.csv到 Pandas 数据框中并查看：

      1. energy = load_data('./data')[['load']]
      2. energy.head(10)

    5. 绘制从 2012 年 1 月到 2014 年 12 月的所有可用能源数据。我们在上一课中看到了这些数据，这应该不足为奇：

      1. energy.plot(y='load', subplots=True, figsize=(15, 8), fontsize=12)
      2. plt.xlabel('timestamp', fontsize=12)
      3. plt.ylabel('load', fontsize=12)
      4. plt.show()

      现在，让我们建立一个模型！

    创建训练和测试数据集

    现在您的数据已加载，因此您可以将其分成训练集和测试集。您将在训练集上训练您的模型。像往常一样，在模型完成训练后，您将使用测试集评估其准确性。您需要确保测试集覆盖训练集中的较晚时间段，以确保模型不会从未来时间段获取信息。

    1. 为训练集分配从 2014 年 9 月 1 日到 10 月 31 日两个月的时间段。测试集将包括 2014 年 11 月 1 日至 12 月 31 日这两个月的时间段：

      1. train_start_dt = '2014-11-01 00:00:00'
      2. test_start_dt = '2014-12-30 00:00:00'

      由于该数据反映了能源的每日消耗量，因此存在强烈的季节性模式，但消耗量与最近几天的消耗量最为相似。

    2. 可视化差异：

      1. energy[(energy.index < test_start_dt) & (energy.index >= train_start_dt)][['load']].rename(columns={'load':'train'}) \
      2.     .join(energy[test_start_dt:][['load']].rename(columns={'load':'test'}), how='outer') \
      3.     .plot(y=['train', 'test'], figsize=(15, 8), fontsize=12)
      4. plt.xlabel('timestamp', fontsize=12)
      5. plt.ylabel('load', fontsize=12)
      6. plt.show()

      因此，使用相对较小的时间窗口来训练数据就足够了。

      注意：由于我们用于拟合 ARIMA 模型的函数在拟合期间使用样本内验证，因此我们将省略验证数据。

    准备训练数据

    现在，您需要通过对数据执行过滤和缩放来准备训练数据。过滤您的数据集以仅包含您需要的时间段和列，并进行缩放以确保将数据投影在 0,1 区间内。

    1. 过滤原始数据集以仅包含上述每组时间段，并且仅包含所需的列“负载”加上日期：

      1. train = energy.copy()[(energy.index >= train_start_dt) & (energy.index < test_start_dt)][['load']]
      2. test = energy.copy()[energy.index >= test_start_dt][['load']]
      3.  
      4. print('Training data shape: ', train.shape)
      5. print('Test data shape: ', test.shape)

      可以看到数据的形状：

      1. Training data shape:  (1416, 1)
      2. Test data shape:  (48, 1)

    2. 将数据缩放到 (0, 1) 范围内。

      1. scaler = MinMaxScaler()
      2. train['load'] = scaler.fit_transform(train)
      3. train.head(10)

    3. 可视化原始数据与缩放数据：

      1. energy[(energy.index >= train_start_dt) & (energy.index < test_start_dt)][['load']].rename(columns={'load':'original load'}).plot.hist(bins=100, fontsize=12)
      2. train.rename(columns={'load':'scaled load'}).plot.hist(bins=100, fontsize=12)
      3. plt.show()

      原始数据

      缩放后的数据

    4. 现在您已经校准了缩放数据，您可以缩放测试数据：

      1. test['load'] = scaler.transform(test)
      2. test.head()

    实施 ARIMA

    是时候实施 ARIMA 了！您现在将使用statsmodels之前安装的库。

    现在您需要执行几个步骤

    1. 通过调用SARIMAX()和传入模型参数来定义模型：p、d 和 q 参数，以及 P、D 和 Q 参数。
    2. 通过调用 fit() 函数为训练数据准备模型。
    3. 调用forecast()函数进行预测并指定要预测的步数 (the horizon)。

    🎓所有这些参数是干什么用的？在 ARIMA 模型中，有 3 个参数用于帮助对时间序列的主要方面进行建模：季节性、趋势和噪声。这些参数是：

    p：与模型的自回归方面相关的参数，它包含过去的值。 d：与模型的积分部分相关的参数，它影响差分量（🎓记住差异👆?) 应用于时间序列。 q：与模型的移动平均部分相关的参数。

    注意：如果您的数据具有季节性方面 - 这个是 - ，我们使用季节性 ARIMA 模型 (SARIMA)。在这种情况下，您需要使用另一组参数：P、D和，它们描述与 、和Q相同的关联p，但对应于模型的季节性分量。dq

    1. 首先设置您的首选水平值。让我们尝试 3 小时：

      1. # Specify the number of steps to forecast ahead
      2. HORIZON = 3
      3. print('Forecasting horizon:', HORIZON, 'hours')

      为 ARIMA 模型的参数选择最佳值可能具有挑战性，因为它有些主观且耗时。您可以考虑使用库auto_arima()中的函数，pyramid

    2. 现在尝试一些手动选择来找到一个好的模型。

      1. order = (4, 1, 0)
      2. seasonal_order = (1, 1, 0, 24)
      3.  
      4. model = SARIMAX(endog=train, order=order, seasonal_order=seasonal_order)
      5. results = model.fit()
      6.  
      7. print(results.summary())

      打印结果表。

    你已经建立了你的第一个模型！现在我们需要找到一种方法来评估它。

    评估您的模型

    要评估您的模型，您可以执行所谓的walk forward验证。在实践中，每次有新数据可用时，都会重新训练时间序列模型。这允许模型在每个时间步进行最佳预测。

    从使用此技术的时间序列开始处开始，在训练数据集上训练模型。然后对下一个时间步进行预测。根据已知值评估预测。然后扩展训练集以包含已知值并重复该过程。

    注意：您应该保持训练集窗口固定以获得更有效的训练，以便每次向训练集中添加新观察值时，都会从集合的开头删除观察值。

    此过程提供了对模型在实践中如何执行的更稳健的估计。然而，它是以创建这么多模型的计算成本为代价的。如果数据很小或模型很简单，这是可以接受的，但在规模上可能是一个问题。

    前向验证是时间序列模型评估的黄金标准，建议用于您自己的项目。

    1. 首先，为每个 HORIZON 步骤创建一个测试数据点。

      1. test_shifted = test.copy()
      2.  
      3. for t in range(1, HORIZON+1):
      4.     test_shifted['load+'+str(t)] = test_shifted['load'].shift(-t, freq='H')
      5.  
      6. test_shifted = test_shifted.dropna(how='any')
      7. test_shifted.head(5)
      loadload+1load+2
      2014-12-3000:00:000.330.290.27
      2014-12-3001:00:000.290.270.27
      2014-12-3002:00:000.270.270.30
      2014-12-3003:00:000.270.300.41
      2014-12-3004:00:000.300.410.57

      数据根据其水平点水平移动。

    2. 使用这种滑动窗口方法在测试数据长度大小的循环中对测试数据进行预测：

      1. %%time
      2. training_window = 720 # dedicate 30 days (720 hours) for training
      3.  
      4. train_ts = train['load']
      5. test_ts = test_shifted
      6.  
      7. history = [x for x in train_ts]
      8. history = history[(-training_window):]
      9.  
      10. predictions = list()
      11.  
      12. order = (2, 1, 0)
      13. seasonal_order = (1, 1, 0, 24)
      14.  
      15. for t in range(test_ts.shape[0]):
      16.     model = SARIMAX(endog=history, order=order, seasonal_order=seasonal_order)
      17.     model_fit = model.fit()
      18.     yhat = model_fit.forecast(steps = HORIZON)
      19.     predictions.append(yhat)
      20.     obs = list(test_ts.iloc[t])
      21.     # move the training window
      22.     history.append(obs[0])
      23.     history.pop(0)
      24.     print(test_ts.index[t])
      25.     print(t+1, ': predicted =', yhat, 'expected =', obs)

      您可以观看正在进行的培训：

      1. 2014-12-30 00:00:00
      2. 1 : predicted = [0.32 0.29 0.28] expected = [0.32945389435989236, 0.2900626678603402, 0.2739480752014323]
      3.  
      4. 2014-12-30 01:00:00
      5. 2 : predicted = [0.3  0.29 0.3 ] expected = [0.2900626678603402, 0.2739480752014323, 0.26812891674127126]
      6.  
      7. 2014-12-30 02:00:00
      8. 3 : predicted = [0.27 0.28 0.32] expected = [0.2739480752014323, 0.26812891674127126, 0.3025962399283795]

    3. 将预测与实际负载进行比较：

      1. eval_df = pd.DataFrame(predictions, columns=['t+'+str(t) for t in range(1, HORIZON+1)])
      2. eval_df['timestamp'] = test.index[0:len(test.index)-HORIZON+1]
      3. eval_df = pd.melt(eval_df, id_vars='timestamp', value_name='prediction', var_name='h')
      4. eval_df['actual'] = np.array(np.transpose(test_ts)).ravel()
      5. eval_df[['prediction', 'actual']] = scaler.inverse_transform(eval_df[['prediction', 'actual']])
      6. eval_df.head()

      输出

      timestamphpredictionactual
      02014-12-3000:00:00t+13,008.743,023.00
      12014-12-3001:00:00t+12,955.532,935.00
      22014-12-3002:00:00t+12,900.172,899.00
      32014-12-3003:00:00t+12,917.692,886.00
      42014-12-3004:00:00t+12,946.992,963.00

      与实际负载相比，观察每小时数据的预测。这有多准确？

    检查模型准确性

    通过测试所有预测的平均绝对百分比误差 (MAPE) 来检查模型的准确性。

    🧮给我看看数学

    MAPE用于将预测精度显示为由上述公式定义的比率。实际t和预测t之间的差异除以实际t。“这个计算中的绝对值是每个预测时间点的总和，然后除以拟合点的数量 n。” 维基百科

    1. 用代码表达方程：

      1. if(HORIZON > 1):
      2.     eval_df['APE'] = (eval_df['prediction'] - eval_df['actual']).abs() / eval_df['actual']
      3.     print(eval_df.groupby('h')['APE'].mean())

    2. 计算一步的 MAPE：

      print('One step forecast MAPE: ', (mape(eval_df[eval_df['h'] == 't+1']['prediction'], eval_df[eval_df['h'] == 't+1']['actual']))*100, '%')

      一步预测 MAPE：0.5570581332313952 %

    3. 打印多步预测 MAPE：

      print('Multi-step forecast MAPE: ', mape(eval_df['prediction'], eval_df['actual'])*100, '%')
      Multi-step forecast MAPE:  1.1460048657704118 %

      一个不错的低数字是最好的：考虑 MAPE 为 10 的预测值下降 10%。

    4. 但与往常一样，这种精度测量更容易直观地看到，所以让我们绘制它：

      1.  if(HORIZON == 1):
      2.     ## Plotting single step forecast
      3.     eval_df.plot(x='timestamp', y=['actual', 'prediction'], style=['r', 'b'], figsize=(15, 8))
      4.  
      5. else:
      6.     ## Plotting multi step forecast
      7.     plot_df = eval_df[(eval_df.h=='t+1')][['timestamp', 'actual']]
      8.     for t in range(1, HORIZON+1):
      9.         plot_df['t+'+str(t)] = eval_df[(eval_df.h=='t+'+str(t))]['prediction'].values
      10.  
      11.     fig = plt.figure(figsize=(15, 8))
      12.     ax = plt.plot(plot_df['timestamp'], plot_df['actual'], color='red', linewidth=4.0)
      13.     ax = fig.add_subplot(111)
      14.     for t in range(1, HORIZON+1):
      15.         x = plot_df['timestamp'][(t-1):]
      16.         y = plot_df['t+'+str(t)][0:len(x)]
      17.         ax.plot(x, y, color='blue', linewidth=4*math.pow(.9,t), alpha=math.pow(0.8,t))
      18.  
      19.     ax.legend(loc='best')
      20.  
      21. plt.xlabel('timestamp', fontsize=12)
      22. plt.ylabel('load', fontsize=12)
      23. plt.show()

    🏆一个非常漂亮的图，显示了一个准确度很高的模型。做得好！

  • 相关阅读:
    农业新闻查询易语言代码
    QT信号槽实现分析
    LeetCode 847. Shortest Path Visiting All Nodes【状态压缩,BFS；动态规划,最短路】2200
    又有新框架上线了，测试、AI 通通有「GitHub 热点速览」
    SpringMVC枚举类型字段处理
    C语言从入门到入土——函数
    超神之路 数据结构 2 —— Queue队列实现和循环队列和普通队列的性能比较
    5. `sklearn`下的线性回归
    Xshell连接显示“服务器发送了一个意外的数据包。received:3,expected:20“问题
    2022.09 青少年Python等级考试（六级） 选择题部分
  • 原文地址：https://blog.csdn.net/sikh_0529/article/details/126961570