时间序列聚类和降维
利用 Kolmogorov Smirnov 统计和机器学习对传感器数据进行聚类
Photo by Igor Ferreira on Unsplash
数据科学家必须小心处理时间序列。这种数据包含关于时间依赖性的内在信息。我们的工作是在可能和有用的地方提取这些黄金资源，以帮助我们的模型发挥最佳性能。
对于时间序列，当我们面对降维或聚类的问题时，我看到了困惑。我们习惯于在更经典的领域中考虑这些任务，而当我们处理时间序列时，它们仍然是一个标签。
在这篇文章中，我试图澄清这些话题，开发一个有趣的解决方案，在这个方案中，我与来自不同个人的多维系列一起工作。我们的目的是利用深度学习，以无监督的方式对它们进行聚类，警惕相关性，并指出每个数据科学家都必须知道的有用技术！
数据集
我从 UCI 机器学习知识库中获得数据；我选择了人体运动基元检测加速度计数据公共数据集。这些数据是标记的加速度计数据记录的公共集合，用于创建和验证人类运动原语的加速度模型。
跟踪不同类型的活动，即喝酒、吃饭、爬山等。对于测量的特定个人的特定活动，我们有 3 个不同的传感器系列可供使用:X 轴(指向手部)、Y 轴(指向左侧)、Z 轴(垂直于手部平面)。
我在这种情况下思考自己，因为它允许在一个单独的案例中执行我们最初的问题聚类(多个个体)和维度缩减(每个个体的多个序列)。
下面我画了两个数据的例子，分别来自男性和女性。总共，我们有 20 个测量长度相同的个体。
降维
首先，我们的攻击计划提供了解决多维度问题。我们希望将存储在传感器数据中的所有信息总结成一个有意义的系列。这一最新步骤将使我们能够轻松地将我们的个人分组。
有很多种技术可以降低数据的维数，但我们的注意力集中在深度学习算法上。一个神经网络结构将允许我们容易地处理我们的初始数据:我记得你，我们有 20 个个体，对于每个个体，我们有 3 个长度为 170 的位置运动序列(用 pythonic 语言来说，我们有一个维数为 20x170x3 的数组)。传统的基于 PCA 的方法不允许我们处理这种问题，所以我们在 Keras 中建立了我们手工制作的自动编码器来处理我们臭名昭著的原始数据结构。
inp = Input(shape=(data.shape[1], data.shape[2]))
encoder = TimeDistributed(Dense(200, activation='tanh'))(inp)
encoder = TimeDistributed(Dense(50, activation='tanh'))(encoder)
latent = TimeDistributed(Dense(10, activation='tanh'))(encoder)
decoder = TimeDistributed(Dense(50, activation='tanh'))(latent)
decoder = TimeDistributed(Dense(200, activation='tanh'))(decoder)out = TimeDistributed(Dense(3))(decoder)autoencoder = Model(inputs=inp, outputs=out)
autoencoder.compile(optimizer='adam', loss='mse')
上面我已经展示了我使用的架构:时间分布层允许处理 3D 数据，其中索引 1 的维度将被认为是时间维度。对于我们的实验，我们使用前 10 个人来训练我们的自动编码器，并利用其余的人来计算相对预测的误差重构。在输入自动编码器之前，不要忘记标准化你的数据！在我们的例子中，我已经通过单个观察(通过行)对每个人的数据进行了标准化。
我们最终的重建误差看起来类似于下面的这个，当我们有信心并且能够检测到这个被选择的人的活动时，我们有接近零的点；而当我们的模型没有足够的学习，并且没有足够的信心来重建行走活动时，我们有很高的价值。
相关聚类
此时，我们已经有了可管理的对象(20 个人的维度为 175x1 ),并且我们已经准备好进行集群。从技术上来说，我们对测试个体的重构误差进行层次聚类。为了捕捉这些系列之间的重要关系，我们尝试了两种不同的工具来组合我们的集群。
在第一阶段，我们的选择包括采用皮尔逊相关指数。不幸的是，这个度量在统计学和机器学习领域非常 过度估计和滥用 ，但是我们想给它一个机会…
获得相关矩阵后，我们直接对其进行操作，执行层次聚类。我们应用高阈值(传感器系列之间的最高成对距离的 99%)来形成我们的扁平聚类。这将导致高水平组的创建，数量很少，但是给我们测试数据的第一印象深刻的概述。
d = sch.distance.pdist(corr)
L = sch.linkage(d, method='ward')
ind = sch.fcluster(L, d.max(), 'distance')
dendrogram = sch.dendrogram(L, no_plot=True)df = [df[i] for i in dendrogram['leaves']]
labels = [person_id[10:][i] for i in dendrogram['leaves']]
corr = np.corrcoef(df)dendrogram = sch.dendrogram(L, labels=[person_id[10:][i] for i in dendrogram['leaves']])
Hierarchical Clustering on Correlation Matrix
查看相关矩阵的颜色强度(我们刚刚对其进行了聚类操作)，我们看不到明显的分组模式。在右边，树状图的切割线(黑线)在一些初始的“不确定性”之后，并没有创造出理性的群体。男女混在一起没有逻辑！
皮尔逊相关指数再次证实了它的不可靠性，我们必须走另一条路…
KOLMOGOROV SMIRNOV 聚类
最近，我读到了 Kolmogorov Smirnov 的统计数据，这对我产生了双重影响:它让我想起了大学，也让我注意到了它的适应性。这个统计量，加上相对 p 值，用于测量两个样本之间的分布差异。我认为我们的聚类任务是这个杀手级工具的一个很好的应用领域。
使用 python 计算这个统计数据非常简单，为了在我们的例子中使用它，我们只需创建 Kolmogorov Smirnov 矩阵(相当于相关矩阵)并重复上面所做的相同步骤。
Hierarchical Clustering on Kolmogorov Smirnov Matrix
现在，查看我们矩阵的颜色强度(我们刚刚对其进行了聚类操作)，我们可以观察到女性和男性之间存在一种模式。从右边的树状图中可以清楚地看到，我们的分级程序已经创建了两个合理的群体，其中男性和女性都是分开的。集群建立之初的“不确定性”也消失了。
这是我们想要的结果，它证实了 Kolmogorov Smirnov 统计数据在每个数据科学家的武器库中的重要性。
摘要
在这篇文章中，我们同时解决了时间序列数据的降维和聚类问题。我们利用自动编码器来总结(以重建误差的形式)加速度计的相关特性。在我们的一维系列中，我们对个体进行了聚类划分。最令人满意的结果来自 Kolmogorov Smirnov 统计和层次聚类的组合，再次证实了必须谨慎处理 Pearson 相关性。
保持联系: Linkedin
时间序列预测(1):初步分析
所有这些信息都可以在我的 kaggle 个人资料中找到。
让我们从定义什么是时间序列开始这个系列。我选择了最简单的方法，并询问了维基百科，答案如下:
时间序列是按时间顺序索引(或列出或绘制)的一系列数据点。最常见的是，时间序列是在连续的等间隔时间点取得的序列。因此，它是一个离散时间数据序列。—维基百科
因此，时间序列基本上与任何其他数据集相似，但具有两个重要特征:
- 数据的顺序很重要，我们不仅关心数据的值，还关心这些值何时出现。
- 我们有观察每个实例的时间信息:我们或者有一个包含日期时间信息的列，或者我们知道数据在时间上是等间距的(例如，每秒一个值)。
在这里，我将介绍时间序列的概念，并对数据进行初步分析和调整；由于 Kaggle 数据通常非常清楚，这项工作比现实生活中容易得多，但其想法是表示过程，而不是编写复杂的代码。我要遵循的流程是这样的:
- 加载数据并设置子集(我不会在这里使用整个数据集)
- 可视化数据
- 清理数据:是否有丢失的值？所有的数据都是期望的格式吗？
- 查看统计数据:有异常值吗？
- 寻找平稳性和自相关性
- 趋势和季节性分解
数据
我正在使用来自 Kaggle 知识库的历史每小时天气数据。该数据集包含不同城市大约 5 年的天气信息，包括温度、湿度、压力、风和天气描述等信息。我将把重点放在温度数据上，因为我觉得它更直观。
为了举例，我将使用两个城市的温度数据，选择的将是旧金山(因为我喜欢它)和变化最大的城市。
加载数据
我将在整个笔记本中使用 Pandas 来处理数据和 matplotlib 进行可视化。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
让我们从加载数据 a 开始，快速浏览一下它的形状和第一个值。因为我已经知道这是一个时间序列，所以我将把 datetime 列设置为索引。
temp = pd.read_csv('../input/temperature.csv', parse_dates=['datetime'])
temp = temp.set_index('datetime')
print('Dataset shape: {}'.format(temp.shape))
temp.head()
该数据集由 36 列组成，对应 36 个不同的城市，45253 行，给出了从 2012 年末到 2017 年末每小时的温度值。
如我所说，我要使用的城市之一是年气温变化最大的城市。让我们找出是哪一个。
all_std = temp.std(axis=0)
max_std = all_std.max()
city_max_std = temp.columns[all_std==max_std][0]print('City with highest temperature variation: {} ({} degrees)'.format(city_max_std,round(max_std,2)))
现在让我们对数据进行子集划分，使其仅包含旧金山和明尼阿波利斯的值，并查看一些统计数据。
data = temp[['San Francisco','Minneapolis']]
data.describe()
我看到的第一件事是两列中都有丢失的数据，我将很快处理这个问题。此外，这个温度值显然不在日常温度的正常范围内。因为我住在欧洲，所以我要将数据从开尔文转换为摄氏度(抱歉，如果你用的是 Farenheit，但是国际系统)。
data = data-273.15
data.describe()
可视化数据
让我们来看看这两个城市历年的气温情况。
_=data.plot(
figsize=(15,5),
subplots=False,
title='Temperature',
alpha=0.7
)
_=plt.xlabel('Date')
_=plt.ylabel('Temperature')
还不错，明尼阿波利斯的天气似乎很冷，也许冬天我不会去那里。
此外，数据有明显的季节性，周期为一年。也有更多的变化，我稍后会检查这种变化是否与一些日常季节性有关，或者它更随机。
从这个图中我可以看出，旧金山的数据更少。蓝线距离 2018 年并没有那么近。
清理数据
正如我们所看到的，这里明显缺少值，图中引起我注意的一些东西可能是原因之一。正如我提到的，旧金山的数据结束得更早，为了同时处理这两个系列，我将丢失明尼阿波利斯系列的最终值。
为此，我将保留旧金山的所有非缺失值，并查看它们到达的最大日期。然后，我们将删除日期大于该日期的所有数据。
SF_non_missing = data['San Francisco'].dropna()
max_date = SF_non_missing.index.max()
data = data[data.index <= max_date]
让我们看看我们是否还有丢失的值。
print(data.isna().sum())
好了，我们还是要处理数据缺失的问题，不过有件事我想先处理一下。我在这里的目的是研究数据的年度行为，我对每日变化不感兴趣，所以我将通过取当天所有温度的平均值、最小值和最大值，将数据重新采样为每日频率。
data_mean = data.resample('D').mean()
data_min = data.resample('D').min()
data_max = data.resample('D').max()
print('Resample shape: {}'.format(data_mean.shape))
data_mean.describe()
如果每天都有更多没有丢失数据的行，这也可以解决我们丢失数据的问题。让我们检查一下。
print('Missing data now?')
print(data_mean.isna().sum())
现在没有丢失数据！这意味着我们每天至少有一个值。
如果不是这样，我会使用前一天的值。对于时间序列，我更喜欢这种解决方案，而不仅仅是删除行(最好在所有行之间有相同的时间间隔)或只使用平均值(因为这会扰乱曲线的形状)。
让我们检查一下重采样后数据的外观。
_=data_mean.plot(
figsize=(15,5),
subplots=False,
title='Temperature',
alpha=0.7
)
_=plt.fill_between(
x=data_mean.index,
y1=data_min['San Francisco'].values,
y2=data_max['San Francisco'].values,
alpha=0.3
)
_=plt.fill_between(
x=data_mean.index,
y1=data_min['Minneapolis'].values,
y2=data_max['Minneapolis'].values,
color='orange',
alpha=0.3
)
_=plt.xlabel('Date')
_=plt.ylabel('Temperature')
曲线周围的阴影显示当天的最小值-最大值，而主线显示平均值。
现在两条曲线都在同一点结束，我们有更少的每日变化。
极端值
有时奇怪的事情会发生，我们最终得到的值会搞乱整个模型。例如，一个传感器可能发生故障，在夏天测量到零下 10 度的温度，这肯定是不正常的。在夏天的其他时候，你可以看到 10，这很低，但可能不是一个错误，有时它会变冷，你不能删除这种类型的值，因为它可能有一些原因。我们必须小心对待异常值，除非我们知道它是一个错误或者是一次性的，不应该影响我们的模型，否则我们不应该删除它。
识别哪些点是异常值以及如何处理它们并不总是很容易，但是一个好的起点是检查数据以查看是否有值非常高或非常低的点。直观地看到这一点的一个好方法是使用直方图。
_=plt.hist(data_mean['San Francisco'], alpha=0.5, label='San Francisco')
_=plt.hist(data_mean['Minneapolis'], alpha=0.5, label='Minneapolis')
_=plt.legend()
让我们分别来看看每个城市:
- 旧金山的值似乎遵循具有小标准偏差的高斯分布，我们在该数据中看不到任何异常值。
- 明尼阿波利斯的曲线不太完美，向右偏高(负偏)。虽然我们不能说这些点中的任何一个是异常值，因为有相当多的异常值，同样，在温度的可视化表示中，我们可以看到每年都达到非常低的值。
我没有看到任何异常值，也不认为我应该删除任何点。
寻找平稳性和自相关性
我们先前绘制的高斯型直方图是时间序列可以是平稳的第一个线索。
另一个线索是计算不同时间范围内的时间序列的一些统计数据，并寻找变化。
cut = data_mean.index[int(0.5*len(data_mean))]
print('Mean before {}:'.format(cut))
print(data_mean.loc[:cut].mean())
print('')
print('Mean after {}:'.format(cut))
print(data_mean.loc[cut:].mean())
print('')
print('---------------------------')
print('')
print('Std before {}:'.format(cut))
print(data_mean.loc[:cut].std())
print('')
print('Std after {}:'.format(cut))
print(data_mean.loc[cut:].std())
我们可以看到，旧金山的值非常接近，但明尼阿波利斯的值更远，它们仍然足够接近，时间序列是平稳的，因为我们需要考虑标准偏差。
这种方法并不证明或否认我们的时间序列是平稳的，它只是表明它可以是平稳的。
我们还可以使用统计测试来确定是否应该拒绝非平稳性假设增强的 Dickey-Fuller 测试。
from statsmodels.tsa.stattools import adfullerresult = adfuller(data_mean['San Francisco'])
print('San Francisco')
print('--------------------------')
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))print('\n\n')
result = adfuller(data_mean['Minneapolis'])
print('Minneapolis')
print('--------------------------')
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
print('\t%s: %.3f' % (key, value))
那么，这些值意味着什么呢？
- ADF 统计量是扩大的 Dicken-Fuller 评分，该值越负，我们可以拒绝零假设(时间序列是平稳的概率)的确定性越高。
- p 值是零假设的置信度。这个值的通常阈值是 0.05，这意味着如果 p 值< = 0.05 ，我们可以拒绝零假设。
- 其余值分别是 99%、95%和 90%置信区间的临界值。
因此，在这种情况下，所有这些意味着我们可以拒绝旧金山的零假设，因为 p 低于 0.05，并且 ADF 得分低于 99%置信区间的限制。然而，我们无法拒绝明尼阿波利斯的这一假设，我们可以说它是具有 90%置信区间的平稳性，但是由于阈值是 95% ( p = 0.05 )，我们不能拒绝它。这意味着我们应该在应用任何模型之前区别数据。
可以在下面找到一些参考资料来理解这些测试:
自相关
最后要检查的是数据是否自相关。我想使用一些自回归方法来预测未来的帖子，我只能在数据自相关的情况下这样做(意味着特定时间点的值取决于以前的值)。
statmodels 库提供了一个很好的工具来检查这一点。阴影区域之外的一切都很有可能是自相关的(超过 95%的置信区间)。
import statsmodels.api as sm
print('San Francisco')
_=sm.graphics.tsa.plot_acf(data_mean['San Francisco'])
plt.show()
print('Minneapolis')
_=sm.graphics.tsa.plot_acf(data_mean['Minneapolis'])
plt.show()
让我们把重点放在最重要的滞后(更接近要点的)，例如，一年范围的数据。
import statsmodels.api as sm
print('San Francisco')
_=sm.graphics.tsa.plot_acf(data_mean['San Francisco'], lags=365)
plt.show()
print('Minneapolis')
_=sm.graphics.tsa.plot_acf(data_mean['Minneapolis'], lags=365)
plt.show()
我们还可以检查部分自相关，它计算去除任何其他先前点(更接近新值的点)的影响的相关性。在这里，更远的点失去了重要性，我将重点放在一个月的范围。
print('San Francisco')
_=sm.graphics.tsa.plot_pacf(data_mean['San Francisco'], lags=30)
plt.show()
print('Minneapolis')
_=sm.graphics.tsa.plot_pacf(data_mean['Minneapolis'], lags=30)
plt.show()
以下是一些关于自相关的参考资料:
趋势-季节性分解
我们可以把时间序列看作趋势、季节性和残差(噪声或其他随机行为)的组合。由这些部分组成的时间序列可以是加法的，也可以是乘法的:
- 加法:数据=趋势+季节性+残差
- 乘法:数据=趋势季节性残差
Statsmodels 包提供了一个函数来一次提取这 3 个组件:季节性分解
这里的分解很简单，因为我们知道周期是 365 天。
from statsmodels.tsa.seasonal import seasonal_decompose as sd
sd_SF = sd(data_mean['San Francisco'], freq=365)
sd_M = sd(data_mean['Minneapolis'], freq=365)_=plt.figure(figsize=(15,10))
ax1=plt.subplot(311)
_=ax1.plot(sd_SF.trend, label='San Francisco', alpha=0.7)
_=ax1.plot(sd_M.trend, label='Minneapolis', alpha=0.7)
_=plt.legend()
ax2=plt.subplot(312)
_=ax2.plot(sd_SF.seasonal, label='San Francisco', alpha=0.7)
_=ax2.plot(sd_M.seasonal, label='Minneapolis', alpha=0.7)
_=plt.legend()
ax3=plt.subplot(313)
_=ax3.plot(sd_SF.resid, label='San Francisco', alpha=0.7)
_=ax3.plot(sd_M.resid, label='Minneapolis', alpha=0.7)
_=plt.legend()
看到这一点，我们可以理解为什么我们发现明尼阿波利斯的数据是不稳定的，温度有明显的上升。
我们还可以通过移动平均来发现趋势，我在这里不会这样做，因为这是季节性分解函数在幕后做的事情。
暂时就这样了。以后的帖子我先说一些预测算法。
希望这能帮助到一些和我一样的初学者，感谢阅读！😃
警察。:我非常希望得到一些反馈
Pd2。:我要感谢 Jason Brownlee 的博客 machinelearningmastery，因为我在这里学到了我所知道的关于时间序列数据的大部分知识，他有许多关于这个主题的帖子，他的解释和例子非常清楚，如果你刚刚开始，或者想扩展你关于这个或任何其他机器学习主题的知识，我完全推荐它
时间序列预测—入门指南
多元时间序列预测导论
介绍
当我开始写这篇文章时，我想解释如何用一个“简单”的时间序列(也就是单变量时间序列)进行预测。然而，我所参与的项目中最具挑战性的部分是预测需要结合多个变量。出于这个原因，我决定让这个指南更接近现实，并使用多元时间序列。
让我们先搞清楚一些概念…
一个 多变量 TS 是一个具有多个时间相关变量的时间序列。每个变量依赖于它过去的值，但也依赖于其他变量。预测值时会考虑这种依赖性。这些变量可以是内生的，也可以是外生的。在这里，我将重点关注外生变量。
外生变量是其值在模型之外确定并强加于模型的变量。换句话说，影响模型而不受其影响的变量。点击阅读更多关于外生变量的信息。
许多模型可以用来解决这样的任务，但我们将与 SARIMAX 合作。SARIMAX 代表带有外生回归量的季节性自回归综合移动平均线。
一切都好！现在，我们将介绍建立销售预测者可以遵循的步骤。
正如我在之前解释的那样，处理时间序列会带来一些挑战，比如让它保持平稳。如果你想知道更多关于我为什么在 dataframe 上执行一些转换的细节，去看看我以前的帖子。本文的重点是预测方法。
在这个机会中，我们有两个文件:一个包含过去销售的数据，另一个包含当地公共假日的信息。可以想象，任务将是通过组合这两个数据集来预测销售额。加载文件后，数据帧看起来像这样:
feriados_df — Holidays dataframe
ventas_df — past sales dataframe
我们的两个数据集的粒度都是在天的级别，也就是说,‘Date’和‘fecha’列都是具有每日频率的索引。如果我们想设置数据集的频率，我们可以运行下面一行:
ventas_df = ventas_df.resample(‘D’).mean() # 'D' for daily frequency
我们将需要连接这两个数据集，以便用我们拥有的所有数据来拟合我们的模型。’ ventas_df '有我们要预测的变量。‘feria dos _ df’包含我们的外生变量。
为了让我们的生活更轻松，最好在加入 feriados_df 之前先站起来 ventas_df。我用来使序列更加平稳的方法包括应用对数变换和差分。固定化系列存储在“ts_log_diff”数据帧中。
test_stationarity(ts_log_diff)
现在我们可以加入 feriados_df 和 ts_log_diff，这是我们改造后的 ventas_df。
data_df = ts_log_diff.join(feriados_df, how='left')
data_df.head()
data_df — joined dataframe
有时在对熊猫进行一些操作后，我们得到的数据帧会丢失频率。为了解决这个问题，我们可以做到:
data_df = data_df.asfreq('D')
是时候来点特色工程了！
人们可以从现有的特性中想出多种创造新特性的方法。为了简单起见，让我们计算下面的列。
-假日 _ 工作日:公共假日是否在工作日
-假日 _ 周末:公共假日是在星期六还是星期天
-是星期几:日期是否在工作日
-是星期几:如果是周末
-在 25 日和 5 日之间:工资通常在这几天发放
**data_df['isweekday']** = [1 if d >= 0 and d <= 4 else 0 for d in data_df.index.dayofweek]
**data_df['isweekend']** = [0 if d >= 0 and d <= 4 else 1 for d in data_df.index.dayofweek]
**data_df['inbetween25and5']** = [1 if d >= 25 or d <= 5 else 0 for d in data_df.index.day]
**data_df['holiday_weekend']** = [1 if (we == 1 and h not in [np.nan]) else 0 for we,h in data_df[['isweekend','Holiday']].values]
**data_df['holiday_weekday']** = [1 if (wd == 1 and h not in [np.nan]) else 0 for wd,h in data_df[['isweekday','Holiday']].values]
让我们对列“Holiday”也应用一次热编码。
data_df = pd.get_dummies(data_df, columns=['Holiday'], prefix=['holiday'], dummy_na=True)
feature-engineered data_df
请我们已经可以预测了！？
亚斯。有点……首先，我们必须将数据分成训练和测试数据。你知道，为了良好的实践和避免过度拟合的东西；)
我们不能仅仅使用 k-folding 方法来将我们的数据集分成训练和测试。这是因为对于 TS，我们必须考虑时间因素。我们可以应用一些技术，其中包括:
- 训练测试分割尊重观察的时间顺序。
- 多重训练测试分割尊重观察的时间顺序。
- 前推验证，每次收到新数据时更新一个模型。
在这种情况下，将使用 1 号。从系列开始到 2019 年 2 月的数据点将用作训练数据。其余的数据点将用于测试。
生成和可视化预测
result_daily = my_train_sarimax(data_df[:'2019-02-28'], i_order=(2,1,2), i_freq='D', i_seasonorder=(2, 1, 1, 12))
在上面的行中，训练数据点’ data _ df[:’ 2019–02–28 ']被传递给函数。值得注意的是，dataframe 中的第一列必须包含要预测的值。其余列是我们的外生变量(即假期和工程特征)。数据帧的频率在“i_freq”参数中给出。参数“i_order”和“i_seasonorder”指定训练模型所需的参数，请查看 SARIMAX 的文档以了解有关这些参数的更多信息。
my_train_sarimax()函数定义如下。
现在是时候验证我们的预测了。为此，我们将使用一个函数来检索预测值，然后将它们与测试数据点中的真实值进行比较。
ypred, ytruth = compare_pred_vs_real(result_daily, data_df, ‘2019–03–01’, exog_validation=data_df[‘2019–03–01’:].iloc[:,1:])
值得一提的是，必须向模型提供要预测的时间范围的外生变量。请记住，这些是模型外部的变量，它需要这些变量来进行预测。
如果我们查看“compare_pred_vs_real()”定义，我们可以看到预测是使用“get_prediction()”函数进行的。可以通过使用“预测平均值”方法来提取这些值。
Performance of our model
我们可以说我们的模型在 MSE 和 RMSE 方面有相当不错的表现。让我们看看我们的预测与实际售出的商品数量相差有多远。
ypred - ytruth
predictions minus true values
但是，等等…为什么我们会看到十进制数值？出售的商品数量必须是整数！
预测的(反)转换
请记住，我们进行了对数转换，然后对数据集应用了差分。为了看到我们的模型估计的实际销售数字，我们必须恢复这些转换。
由于差异操作，TS 中的原始第一个日期丢失。我们需要从“data_df”中填充缺失的值。接下来，我们需要将预测之前的所有日期追加到‘y _ pred’中。这些日期也来自‘data _ df’。完成所有这些后，我们可以用 cumsum()恢复差异，然后应用 exp()恢复对数转换。
#create a series with the dates that were dropped with differencing
restore_first_values = pd.Series([6.008813], index=[pd.to_datetime(‘2014–01–01’)])#get the values that the prediction does not have
missing_part = data_df[‘cantidad’][:’2019–02–28']
rebuilt = restore_first_values.append(missing_part).append(ypred)#revert differencing:
rebuilt = rebuilt.cumsum()#revert log transformation:
rebuilt = np.exp(rebuilt).round() # apply round() to have integers
我们终于可以看到我们的预测值，并与实际值进行比较。耶！
# Check how far were the predictions from the actual values
rebuilt['2019-03-01':] - ventas_df['cantidad']['2019-03-01':]
Looks like we got it right most of the time 😃
最终意见
请记住，可以使用许多方法来实现 TS 中的平稳性。此外，SARIMAX 不是唯一一个对时间序列进行预测的模型，进一步的参数调整有助于提高模型的准确性。
请随意在我的 github repo 上查看该指南的全部代码。
感谢阅读！
英国道路交通事故的时间序列预测
在这篇文章中，我们旨在预测英国未来的交通事故数量。使用的数据来自交通部(GB) [1]。该数据提供了 2014 年至 2017 年英国人身伤害道路事故情况的详细道路安全信息。我们将使用 ARIMA 和 Prophet 实现时间序列预测，以找出未来道路交通事故的数量。
Photo by Stephen Dawson on Unsplash
什么是时间序列？
时间序列是变量在不同时间取值的一组观察值。例如销售趋势、股票市场价格、天气预报等。时间序列是用来根据以前的观察值预测未来值的。
时间序列的成分
- 趋势:趋势可以显示一个时间序列在很长一段时间内的增长或下降。这种趋势将持续很长一段时间。例如，价格和进出口数据反映出明显的增长趋势。
- 季节性:这些是由于季节性因素而在数据中出现的短期变动。短期通常被认为是随着天气或节日的变化，时间序列发生变化的时期
- 不规则性:这些是在时间序列中发生的不太可能重复的突然变化。它们是时间序列的组成部分，无法用趋势、季节或周期运动来解释。这些变化有时被称为残差或随机分量。
- 周期性的:这些是发生在时间序列中的长期振荡。这些振荡大多在经济数据中观察到，并且这种振荡的周期通常从五年延长到十二年或更长。这些波动与众所周知的商业周期有关。[2]
什么是 ARIMA 模式？
ARIMA 代表自回归综合移动平均线。有季节性和非季节性 ARIMA 模型可用于预测。ARIMA 模型的特征在于 3 项:p，d，q，其中 p 是 AR 项的阶，q 是 MA 项的阶，d 是使时间序列稳定所需的差的数量。如果一个时间序列有季节模式，那么你需要添加季节项，它就变成了萨里玛，是“季节性 ARIMA”的缩写。等我们结束 ARIMA 后会有更多的报道。[3]
让我们用所需的库准备好我们的环境，然后导入数据！
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
import matplotlib
plt.style.use('ggplot')
import warnings
import itertools
matplotlib.rcParams['axes.labelsize'] = 14
matplotlib.rcParams['xtick.labelsize'] = 12
matplotlib.rcParams['ytick.labelsize'] = 12
matplotlib.rcParams['text.color'] = 'k'
检查数据
df = pd.read_csv('~/accident_UK.csv')
df.head()
df.info()
将日期列转换为日期类型
df['Date'] = pd.to_datetime(df['Date'])
df.head()
按日期对数据进行排序
df = df.sort_values(by=['Date'])
df.head()
设置索引的日期
accident = df.set_index('Date')
accident.index
我们来提取一下每个月的平均事故次数。
y = accident['Total_Accident'].resample('MS').mean()
y.head()
让我们看一下每个月平均事故的数量。
y.plot(figsize=(15, 6))
plt.show()
让我们使用时间序列分解来可视化数据，这允许我们将时间序列分解为三个不同的部分:趋势、季节性和噪声。
from pylab import rcParams
import statsmodels.api as sm
rcParams['figure.figsize'] = 16, 10
decomposition = sm.tsa.seasonal_decompose(y, model='additive')
fig = decomposition.plot()
plt.show()
季节性 ARIMA 的参数组合示例
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
print('Examples of parameter combinations for Seasonal ARIMA...')
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[1]))
print('SARIMAX: {} x {}'.format(pdq[1], seasonal_pdq[2]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[3]))
print('SARIMAX: {} x {}'.format(pdq[2], seasonal_pdq[4]))
这一步是为我们的 ARIMA 时间序列模型选择参数。我们在这里的目标是使用“网格搜索”来找到为我们的模型产生最佳性能的最佳参数集。
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(y,
order=param,
seasonal_order=param_seasonal,
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
except:
continue
上述输出表明，SARIMAX (1，1，1)x(1，1，0，12)产生的 AIC 值最低，为 186.26。因此，我们应该认为这是一个最佳选择。
拟合 ARIMA 模型
mod = sm.tsa.statespace.SARIMAX(y,
order=(1, 1, 1),
seasonal_order=(1, 1, 0, 12),
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
print(results.summary().tables[1])
我们需要始终运行模型诊断来调查任何不寻常的行为。
results.plot_diagnostics(figsize=(16, 8))
plt.show()
预测评估
为了了解我们预测的准确性，我们将预测的事故数量与时间序列的实际事故数量进行比较，我们将预测从 2017 年 1 月 1 日开始到数据结束。
pred = results.get_prediction(start=pd.to_datetime('2017-01-01'), dynamic=False)
pred_ci = pred.conf_int()
ax = y['2014':].plot(label='observed')
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7, figsize=(14, 7))
ax.fill_between(pred_ci.index,
pred_ci.iloc[:, 0],
pred_ci.iloc[:, 1], color='k', alpha=.2)
ax.set_xlabel('Date')
ax.set_ylabel('Furniture Sales')
plt.legend()
plt.show()
让我们找出 MSE 来看看我们模型的准确性。均方误差(MSE)主要用作确定算法性能的度量。此外，MSE 是一个变量的观察值和预测值之差的平方的平均值。
y_forecasted = pred.predicted_mean
y_truth = y['2017-01-01':]
mse = ((y_forecasted - y_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))
可视化预测
正如我们在下面的图表中所看到的，英国的交通事故数量在未来几年将会下降。
pred_uc = results.get_forecast(steps=100)
pred_ci = pred_uc.conf_int()
ax = y.plot(label='observed', figsize=(14, 7))
pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
ax.fill_between(pred_ci.index,
pred_ci.iloc[:, 0],
pred_ci.iloc[:, 1], color='k', alpha=.25)
ax.set_xlabel('Date')
ax.set_ylabel('Furniture Sales')
plt.legend()
plt.show()
什么是先知模型？
Prophet 是一个来自脸书的开源时间序列预测算法，它的设计是为了在没有时间序列预测或统计方面的专家知识的情况下易于使用。时间序列预测通过寻找最佳平滑线来建立模型，该最佳平滑线可以表示为以下分量的总和:
- 总体增长趋势
- 早期季节性
- 每周季节性
- 节日影响
先知方法的好处:
- 数据之间不均匀的时间间隔不是问题
- 天娜不是问题
- 默认情况下，会处理多个期间(周和年)的季节性
- 默认设置下工作良好，参数易于解释
让我们按日期对值进行排序
df = df.sort_values(by=['Date'])
df.head()
Prophet 要求时间序列中的变量名为:
- y-目标
- ds —日期时间
因此，下一步是根据上述规范转换数据帧
df = df.rename(columns={'Date': 'ds',
'Total_Accident': 'y'})
df.head()
让我们想象一下每天的交通事故数量
ax = df.set_index('ds').plot(figsize=(15, 8))
ax.set_ylabel('Total Accident')
ax.set_xlabel('Date')
plt.show()
拟合先知模型
将不确定性区间设置为 95%(Prophet 默认为 80%)
from fbprophet import Prophet
my_model = Prophet(interval_width=0.95)
my_model.fit(df)
要使用我们的模型创建预测，我们需要创建一些未来日期。Prophet 为我们提供了一个辅助函数，叫做 make_future_dataframe。我们传入未来周期的数量和频率。以上是我们对未来 36 个月或 3 年的预测。
future_dates = my_model.make_future_dataframe(periods=36, freq='MS')
future_dates.tail()
正如我们在下表中看到的，yhat 是我们的预测值。
forecast = my_model.predict(future_dates)
forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()
让我们创建一个实际值和预测值的图表
plt.figure(figsize=(10,8))
my_model.plot(forecast,
uncertainty=True)
plot_components 为我们提供了趋势和季节性的图表
my_model.plot_components(forecast)
模型评估
from fbprophet.diagnostics import cross_validation
df_cv = cross_validation(my_model, initial='730 days', period='180 days', horizon = '365 days')
df_cv.head()
from fbprophet.diagnostics import performance_metrics
df_p = performance_metrics(df_cv)
df_p.head()
非常感谢你阅读我的文章！！！
参考资料:
[1].https://data . gov . uk/dataset/CB 7 AE 6 f 0-4be 6-4935-9277-47 E5 ce 24 a 11 f/road-safety-data
[2].克利夫兰，R. B .，克利夫兰，W. S .，麦克雷，J. E .，&特彭宁，国际法院(1990)。STL:基于黄土的季节趋势分解过程。官方统计杂志， 6 (1)，3–33。
[3].Hyndman，R. J .，& Khandakar，Y. (2008 年)。自动时间序列预测:R. 统计软件学报， 27 (1)，1–22 的预测包。
时间序列预测:从天真到 ARIMA 及以后
简单和复杂的商业决策预测技术
Photo by Aron Visuals on Unsplash
- 在未来五年内再建一座发电厂需要预测未来的需求；
- 下周在呼叫中心安排员工需要预测下周的呼叫量；
- 储存零售存货需要预测存货需求。
事实上，短期和中期预测是各行各业商业决策的重要组成部分。历史数据是这一预测过程的重要输入。
时间序列数据集是在任何业务中最广泛生成和使用的数据。它们既用于理解过去，也用于预测未来。在本文中，我将通过相关的案例研究和示例来讨论如何将各种预测方法应用于时间序列数据集。这篇文章的目的不是评价哪个模型是好是坏，而是展示我们在实践中可以做预测的许多不同的方法。本文分为两大部分:
(1)首先，我将概述时间序列数据以及如何分解不同的时间序列成分；
(2)然后我将提供不同预测技术的例子以及相关的实现方法
理解和分解时间序列数据
那么什么是时间序列数据呢？时间序列是在离散时间点记录的一系列观察值。它们可以以每小时(如气温)、每天(如 DJI 平均)、每月(如销售额)或每年(如国内生产总值)的时间间隔记录。时间序列广泛应用于所有学科，包括统计学、数学、天文学和工程学。时间序列绘图是在数据科学和分析的探索性数据分析(EDA)阶段实施的最基本的绘图练习之一。
Time series data: Daily Dow Jones Industrial (DJI) averages over 6 months
我们来看看下图。
即使用肉眼也很少有明显的特性:
- 这是一个时间序列图
- 有上升的趋势
- 观察值具有季节性(在固定的时间间隔内起伏)
事实上，如果我们分解这些数据，我们会看到这些组件，如下所示。
时间序列数据的预测技术
预测未来可以简单到外推历史观察的趋势，应用复杂的算法。在这里，我介绍了商业应用程序中最常用的不同技术。我使用两个统计和可视化包在 R 环境中演示这些技术:
# install packages
library("forecast")
library("ggplot2")
在导入数据(零售销售数据集)并将其转换为时间序列对象后，它看起来是这样的。
# import dataset
retaildata = readxl::read_excel("../retail.xlsx", skip=1)# convert dataset into a time series object
retailts = ts(retaildata, start = c(1982,4), frequency = 12 )# plot the time series object with autoplot() function that comes with forecast() package
autoplot(retailts)
首先我将解释一些简单的技术来预测这个时间序列对象。虽然这些简单的技术在实践中不经常使用，但知道它们是如何产生的是很好的，然后我们将进入更复杂的技术。
简单的预测技术
简单的预测技术被用作基准。它们提供了对历史数据的一般理解，并建立直觉，在此基础上增加额外的层复杂性。几种这样的技术在文献中很常见，如:均值模型、简单预测、随机游走、漂移法等。例如，均值模型采用以前观察值的均值，并将其用于预测。那就是:
预测值=训练数据的平均值
另一方面，随机游走将预测下一个值，ŷ(t，它等于前一个值加上一个常数变化。
ŷ(t)= y(t-1)+α
初始预测的值是基于最后一次观察的值设置的。和其他简单的方法一样，它提供了一个大概的数字作为初步估计，直到进一步的研究完成。下面是季节性朴素预测模型的实现。
# subsetting train and test data using window()function
retailtstrain = window(retailts, end = c(2010,3))
retailtstest = window(retailts, start = c(2010,4))# model building
retailtstrain_snaive = snaive(retailtstrain, h=50)# ploting
autoplot(retailts)+
autolayer(retailtstrain_snaive, series = "S Naive", PI = FALSE)
# Accuracy test
accuracy(retailtstrain_snaive, retailtstest)
指数平滑法
根据数据类型的不同，指数平滑几乎没有变化。简单的指数平滑法用于没有明确趋势的非季节性数据，而霍尔特-温特法用于有趋势和季节性的数据。它应用于平稳时间序列，其中平滑由参数 alpha (0~1)控制，其中较低的值意味着对最近观察值的较低权重。下面是使用相同数据集的霍尔特-温特指数平滑的实现。
# Holt Winter Exponential Smoothing
# Two variations: additive for roughly constant seasonal variation, otherwise multiplicative methoddata = window(retailts, start=2000)
data_add = hw(data, seasonal="additive")
data_mult = hw(data, seasonal="multiplicative")autoplot(data, series = "original data")+
autolayer(data_add, series = "additive", PI=FALSE)+
autolayer(data_mult, series = "multiplicative", PI=FALSE)
# model summary and performance
data_add[["model"]]
ARIMA 家族
自回归综合移动平均(ARIMA)可以说是最流行和最广泛使用的统计预测技术。顾名思义，这个技术家族有 3 个组成部分:a)一个“自回归”部分，它模拟序列和它的滞后观测值之间的关系；b)将预测建模为滞后预测误差的函数的“移动平均”模型；和 c)使系列静止的“集成”组件。
该模型采用以下参数值:
定义滞后数量的 p；
d 指定使用的差异数；和
q 定义了移动平均窗口的大小
使用时序对象进行预测的 ARIMA 模型的实现如下:
# implementing auto.arima() to forecast
retailts_arima = auto.arima(retailts, seasonal=TRUE, stepwise = FALSE, approximation = FALSE)
retailts_arima %>% forecast(h=10) %>% autoplot()
更先进的技术
最后，还有一些盒子外的技术，如基于代理的和系统动态建模。系统动力学是 20 世纪 50 年代由麻省理工学院斯隆管理学院开发的，是一种模拟复杂系统行为的方法，其中一个组件的变化会导致其他组件的变化。这种方法广泛应用于医疗保健、疾病研究、公共交通、商业管理和收入预测等行业。系统动力学最著名的应用可能是罗马俱乐部的增长极限模型。
一个系统动态模型代表了一个复杂的系统，通过反馈循环来预测系统的行为。假设一个银行账户有 100 美元的“股票”。每月存入 20 美元(由“流程”1 表示)，每月提取 15 美元(由“流程”2 表示):
System dynamic modeling workflow
在这个例子中，股票的未来价值(即账户存款)被建模为未来流量(即存款和取款)的函数。本文提供了在建筑行业实施系统动态预测模型的案例研究。
额外资源
- 罗布·海德曼的网站https://robjhyndman.com/是 R 中所有预测问题的一站式解决方案。也可以查看开放存取书籍https://otexts.com/fpp2/中的理论讨论以及 R 实施示例。
- 《R 对于时间序列的一个小本子》什么都有简化形式。http://www . Calvin . edu/~ stob/courses/m344/S15/a-little-book-of-r-for-time-series . pdf
- 如果你喜欢 Python，请点击下面的链接。还有一本关于时序的 python 实现的书的参考。https://machinelementmastery . com/time-series-forecasting-methods-in-python-cheat-sheet/
基于深度堆叠单向和双向 LSTMs 的时间序列预测
这篇文章假设读者对 LSTMs 的工作原理有一个基本的了解。不过，你可以在这里得到 LSTMs 的简要介绍。另外，如果你是时间序列预测的绝对初学者，我推荐你去看看这个博客。
这篇文章的主要目的是展示深度堆叠单向和双向 LSTMs 如何作为基于 Seq-2-Seq 的编码器-解码器模型应用于时间序列数据。我们将首先介绍这个架构，然后看看实现它的代码。
让我们深入研究这个模型
等等。！在进入架构的细节之前，让我们理解序列到序列模型的特殊性。因此，顾名思义，序列对序列模型将一系列特征作为输入，并输出一个目标序列作为输入目标序列的延续(它可以预测未来的“n”个时间步)。
它基本上有两个部分，编码器输出输入序列的上下文向量(编码)，然后传递给解码器解码和预测目标。
Lol，这可能有点让人不知所措，但是随着我们进一步深入和可视化架构，你会慢慢理解这些术语的。
模型架构
让我们从基本的编码器-解码器架构开始，然后我们可以逐步添加新的功能和层来构建更复杂的架构。
1。使用单向 LSTMs 作为编码器
这里，LSTM 编码器将时间序列作为输入(每个 LSTM 单元一个时间步长)，并创建输入序列的编码。该编码是由所有编码器 LSTM 单元的隐藏和单元状态组成的向量。编码然后作为初始状态与其他解码器输入一起被传递到 LSTM 解码器，以产生我们的预测(解码器输出)。在模型训练过程中，我们将目标输出序列设置为模型训练的解码器输出。
2。使用双向 LSTMs 作为编码器
双向 LSTMs 具有两个递归分量，前向递归分量和后向递归分量。前向组件计算隐藏和单元状态，类似于标准单向 LSTM，而后向组件通过以逆时间顺序(即从时间步长 Tx 到 1 开始)取输入序列来计算它们。使用后向组件的直觉是，我们正在创建一种方式，网络可以看到未来的数据，并相应地学习其权重。这可能有助于网络捕获标准(前向)LSTM 无法捕获的一些相关性。BLSTM 也是大多数 NLP 任务的首选算法，因为它能够很好地捕捉输入序列中的相关性。
在 BLSTMs 中，前向组件的隐藏和单元状态与后向组件的不同。因此，为了获得编码，前向组件的隐藏和单元状态必须分别与后向组件的隐藏和单元状态连接。
3。使用堆叠单向 LSTMs 作为编码器
当这些层堆叠在一起时，编码器和解码器的第一层 LSTM 单元的输出(单元状态)被传递到第二层 LSTM 单元作为输入。似乎具有几个隐藏层的深度 LSTM 架构可以有效地学习复杂的模式，并且可以逐步建立输入序列数据的更高级别表示。
双向 LSTMs 也可以以类似的方式堆叠。第一层的前向和后向分量的输出分别传递给第二层的前向和后向分量。
实施细节:D
数据准备—
我们将在标准的“墨尔本每日最低气温”(单变量时间序列)数据集上应用上述模型(从这里下载)。但是，在进入培训场景之前，让我们首先准备数据。
- 将数据标准化。
注意: reshape 仅用于将单变量 1D 数组转换为 2D，如果数据已经是 2D，则不需要调用。
2.初始化参数。
3.生成输入输出序列对。
上述函数返回一批大小为“total_start_points”的输入序列和输出(目标)序列，它们将分别被馈送到编码器和解码器。注意返回的序列是形状的 3D 张量(batch_size，input_seq_len，n_in_features ),因为 keras 要求输入为这种格式。
模型—
我们定义了一个单一的函数来构建架构，这取决于传递的隐藏维度列表，以及在调用该函数时设置为“真”或“假”的参数“双向”。
注意当编码器是双向的时，我们在解码器 LSTM 中有‘hidden _ dim * 2’来容纳级联的编码器状态。
模特.飞度:D
默认情况下，Keras 在训练时会打乱数据，因此我们可以(不一定)在“model.fit”函数中设置“shuffle=False ”,因为我们已经在随机生成序列。
注意我们输入零作为解码器输入，也可以使用 教师强制 (其中一个解码器单元的输出作为输入输入到下一个解码器单元)(此处未涉及)。
结果
Unidirectional Layer-1 and Layer-2
Bidirectional Layer-1 and Layer-2
哇！所有的模型只被训练了 100 个时期(具有相同的参数),与单向 lstm 相比，双向 lstm 在学习数据中的复杂模式方面表现突出。因此，所描述的模型可以应用于许多其他时间序列预测方案，甚至可以应用于多变量输入情况，在这种情况下，您可以将具有多个特征的数据作为 3D 张量进行传递。
你可以在我的 GitHub 库中找到这个例子的 Jupyter 笔记本实现。
我希望你喜欢这篇文章，并让你很好地理解了如何使用深度堆叠 LSTMs 进行时间序列预测。非常感谢您的反馈或改进建议。
使用 TensorFlow 2 和 Keras 在 Python 中使用 LSTMs 进行时间序列预测
利用 LSTMs 进行时间序列预测的数据准备和预测介绍
TL；DR 了解时间序列，并使用递归神经网络进行预测。准备序列数据并使用 LSTMs 进行简单预测。
想学习如何使用多元时间序列数据？阅读下一部分:
[## 使用 TensorFlow 2 和 Keras 在 Python 中使用 LSTMs 进行需求预测
了解如何通过深度学习从多元时间序列数据中预测需求
towardsdatascience.com](/demand-prediction-with-lstms-using-tensorflow-2-and-keras-in-python-1d1076fc89a0)
通常，您可能不得不处理包含时间成分的数据。不管你怎么眯眼，都很难做出你喜欢的数据独立性假设。似乎数据中的新值可能依赖于历史值。你怎么能使用这种数据来建立模型呢？
本指南将帮助您更好地理解时间序列数据，以及如何使用深度学习(递归神经网络)建立模型。您将学习如何预处理时间序列，构建一个简单的 LSTM 模型，训练它，并使用它来进行预测。以下是步骤:
- 理解什么是时间序列
- 了解递归神经网络
- 在 Keras 中用 LSTMs 预测时间序列数据
- 评估模型
时间序列
时间序列是数据点的集合，根据它们被收集的时间进行索引。大多数情况下，数据是以固定的时间间隔记录的。时间序列数据的特殊之处是什么？
预测未来时间序列值是实践中一个相当常见的问题。预测下一周的天气、明天的比特币价格、圣诞节期间你的销售数量以及未来的心脏衰竭都是常见的例子。
时间序列数据引入了对先前时间步长的“硬依赖”,因此观测值独立性的假设不成立。时间序列可以具有哪些属性？
平稳性、季节性和自相关性是您可能感兴趣的时间序列的一些属性。
当均值和方差随时间保持不变时，称时间序列为平稳。如果平均值随时间变化，那么时间序列具有趋势。通常，您可以通过应用对数变换来消除它并使序列平稳。
季节性是指特定时间范围内的变化现象。人们在圣诞节期间购买更多的圣诞树(谁会想到)。消除季节性的一种常见方法是使用差异。
自相关 指当前值与前一时间(滞后)的拷贝值之间的相关性。
为什么我们想要季节性、趋势性和平稳的时间序列？这是用经典方法如 ARIMA 模型进行时间序列预测所需的数据预处理步骤。幸运的是，我们将使用递归神经网络进行建模。
递归神经网络
递归神经网络(RNNs)可以预测序列中的下一个值或对其进行分类。一个序列被存储为一个矩阵，其中每一行都是一个描述它的特征向量。自然，矩阵中行的顺序很重要。
rnn 非常适合解决自然语言处理(NLP)任务，其中文本中的单词形成序列，它们的位置很重要。也就是说，前沿 NLP 使用变压器来完成大多数(如果不是全部)任务。
你可能已经猜到了，时间序列只是序列的一种类型。我们必须将时间序列切割成更小的序列，这样我们的 RNN 模型就可以用它们进行训练。但是我们如何训练 rnn 呢？
首先，让我们对递归的含义有一个直观的理解。rnn 包含循环。每个单元都有一个状态，并接收两个输入-来自前一层的状态和来自前一时间步的该层的统计数据。
反向传播算法在应用于 RNNs 时会因为循环连接而失效。展开网络可以解决这个问题，在网络中，具有循环连接的神经元的副本被创建。这将 RNN 转换成常规的前馈神经网络，并且可以应用经典的反向传播。该修改被称为通过时间的反向传播。
经典 rnn 的问题
展开的神经网络会变得非常深(他就是这么说的)，这给梯度计算带来了问题。权重可以变得非常小(消失梯度问题)或者非常大(爆炸梯度问题)。
经典的 rnn 也有记忆问题(长期依赖性)。由于最近状态的压倒性影响，我们用于训练的序列的乞求倾向于被“遗忘”。
在实践中，这些问题可以通过使用门控 rnn 来解决。它们可以存储信息以备后用，就像有一个内存一样。阅读、写作和从记忆中删除都是从数据中学习的。两种最常用的门控 rnn 是长短期记忆网络和门控复发性单位神经网络。
用 LSTMs 进行时间序列预测
我们将从一个简单的例子开始，使用简单的 LSTM 网络预测正弦函数的值。
设置
让我们从库导入和设置种子开始:
数据
我们将从正弦函数中生成 1，000 值，并将其作为训练数据。但是，我们会添加一点精加工:
从正态分布中提取的随机值被添加到每个数据点。那会让我们模型的工作变得有点困难。
数据预处理
对于我们的模型，我们需要将数据“切碎”成更小的序列。但首先，我们将把它分为训练和测试数据:
800 200
为时间序列预测(尤其是 LSTMs)准备数据可能很困难。直观地说，我们需要利用历史( n 时间步长)来预测当前时间步长的值。这里有一个通用函数可以完成这项工作:
该函数的妙处在于它可以处理单变量(单特征)和多变量(多特征)时间序列数据。让我们用 10 个时间步骤的历史来制作我们的序列:
(790, 10, 1) (790,)
我们有形状为(samples, time_steps, features)的序列。我们怎么用它们来做预测呢？
建模
在喀拉斯培养 LSTM 模式很容易。我们将使用序列模型中的 LSTM 层进行预测:
LSTM 图层期望时间步长数和要素数能够正常工作。模型的其余部分看起来像一个常规的回归模型。我们如何训练一个 LSTM 模特？
培养
训练时间序列模型时要记住的最重要的事情是不要打乱数据(数据的顺序很重要)。其余的都很标准:
我们的数据集非常简单，包含了我们采样的随机性。经过大约 15 个时期后，模型已经基本完成了学习。
估价
让我们从我们的模型中选取一些预测:
我们可以根据时间序列的真实值绘制预测图:
我们的预测在这个尺度上看起来非常好。让我们放大:
该模型似乎在捕捉数据的一般模式方面做得很好。它未能捕捉到随机波动，这是一件好事(它概括得很好)。
结论
恭喜你！你做了你的第一个递归神经网络模型！您还了解了如何预处理时间序列数据，这是一件让很多人感到困惑的事情。
我们只是触及了时间序列数据和如何使用递归神经网络的表面。一些有趣的应用是时间序列预测、(序列)分类和异常检测。有趣的部分才刚刚开始！
想学习如何使用多元时间序列数据？阅读下一部分:
[## 使用 TensorFlow 2 和 Keras 在 Python 中使用 LSTMs 进行需求预测
了解如何通过深度学习从多元时间序列数据中预测需求
towardsdatascience.com](/demand-prediction-with-lstms-using-tensorflow-2-and-keras-in-python-1d1076fc89a0)
参考
Scikit-Learn、TensorFlow 和 Keras 深度学习实践指南了解如何解决现实世界的机器学习…
leanpub.com](https://leanpub.com/Hackers-Guide-to-Machine-Learning-with-Python)
原载于https://www.curiousily.com。
使用 LSTM 预测序列数据:简介
使用长短期记忆(LSTM)网络预测股票等序列数据的未来。
Photo by Chris Liverani on Unsplash
预测是使用当前和以前的数据预测未来的过程。主要的挑战是理解数据序列中的模式，然后使用这种模式来分析未来。如果我们手工编码这些模式，那么对于下一个数据来说，这将是乏味的和变化的。深度学习已被证明在理解结构化和非结构化数据中的模式方面更好。
为了理解一长串数据中的模式，我们需要网络来分析时间上的模式。递归网络通常用于学习这种数据。他们能够理解长期和短期的依赖性或时间差异。
好了，不介绍了…
这篇文章将向你展示如何在 Keras 中使用 LSTM 网络和一些很酷的可视化工具来实现一个预测模型。我们将使用来自雅虎财经的谷歌股价，但你可以随意使用任何你喜欢的股票数据。
履行
我已经使用 colab 来实现这个代码，以使可视化更容易，你可以使用你喜欢的方法。我们将从导入必要的库开始:
下载完您的。csv 文件，使用 pandas 加载数据。
dataframe 的信息显示如下:
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 3797 entries, 0 to 3796
Data columns (total 7 columns):
Date 3797 non-null object
Open 3797 non-null float64
High 3797 non-null float64
Low 3797 non-null float64
Close 3797 non-null float64
Adj Close 3797 non-null float64
Volume 3797 non-null int64
dtypes: float64(5), int64(1), object(1)
memory usage: 207.7+ KB
对于本教程，我们只需要 Date 和 Close 列，其他的都可以删除。
在我们进行训练和预测之前，让我们看看数据是什么样子的。对于所有的可视化，我使用的是 Plotly python 库。为什么是 Plotly？…因为它是最好的图形库，可以生成一些好看的图形。
使用 plotly，我们可以定义一个轨迹和布局，它可以做任何其他事情。
图表从 2018 年开始振荡，序列并不平滑…继续前进。
数据预处理
对于我们的分析，让我们在前 80%的数据上训练模型，并在剩余的 20%上测试它。
在我们进行培训之前，我们需要对我们的数据进行一些重大修改。记住，我们的数据仍然是一个序列…一列数字。神经网络被训练为监督模型。因此，我们需要将序列数据转换成监督数据😱。
我来解释一下，训练任何机器学习模型的神经网络都需要数据是{ <特征 >，< 目标 > }格式。同样，我们需要将给定的数据转换成这种格式。这里，我们引入一个回望的概念。
回顾只不过是使用前几天的数据来预测第二天的值。比如我们说回望是 2；所以为了预测明天的股价，我们需要今天和昨天的股价。
回到格式，在给定的一天x(t),特征是 x(t-1)，x(t-2)，…。，x(t-n) 其中 n 为回看。
所以如果我们的数据是这样的，
[2, 3, 4, 5, 4, 6, 7, 6, 8, 9]
所需的数据格式(n=3)如下:
[2, 3, 4] -> [5]
[3, 4, 5] -> [4]
[4, 5, 4] -> [6]
[5, 4, 6] -> [7]
[4, 6, 7] -> [6]
[6, 7, 6] -> [8]
[7, 6, 8] -> [9]
咳😓。
幸运的是，Keras 中有一个模块可以做到这一点: TimeseriesGenerator 。请查阅文档以获取更多信息。
我已经将 look_back 设置为 15，但是您可以随意使用这个值。
神经网络
现在我们的数据已经准备好了，我们可以继续创建和训练我们的网络。
使用 Adam 优化器和均方损失函数训练 25 个时期的 LSTM 单元的简单架构。注意，我们没有使用 model.fit()，而是使用 model.fit_generator()，因为我们已经创建了一个数据生成器。
想了解更多关于 LSTM 网络的信息，请看这篇令人敬畏的博文。
预言；预测；预告
现在我们已经完成了培训，让我们看看网络是否运行良好。我们可以在测试数据上测试模型，看看预测值和实际值是否重叠。
我们可以画出预测值和实际值之间的损失，而不是计算它们。
从图表中，我们可以看到预测值和实际值(地面实况)有些重叠。但如果放大看，拟合并不完美。我们应该预料到这一点，因为这是不可避免的，因为我们正在进行预测。
预测
我们的测试表明这个模型还不错。所以我们可以继续预测未来。
伏笔:因为我们试图预测未来，所以预测中会有很大的不确定性。
预测未来很容易…要预测明天的价值，将过去 n (回望)天的价值输入模型，我们就可以得到明天的价值作为输出。要获得后天的值，需要输入过去 n-1 天的值，以及明天的值和后天的模型输出。
更长时间的预测是不可行的。所以，让我们预测一个月的股票价格。
现在描绘未来的价值，
预测未来很容易…是吗？
在预测未来时，很有可能模型输出在很大程度上是不确定的。模型的输出作为输入反馈给它。这导致模型的噪音和不确定性被重复和放大。
但是，我们还是创建了一个模型，它给出了图表的趋势，以及未来可能的取值范围。
结论
通过这个模型，我们创建了一个能够在一定程度上进行预测的基本模型。尽管这个模型并不完美，但我们有一个能很好地逼近过去数据的模型。但是对于新数据，我们需要更多的参数调整。
如果在模型训练期间存在任何问题，请尝试调整这些参数:
- 回望 _ 回顾
- 批量
- LSTM 台
- 次数
您还可以通过堆叠 LSTM 图层来增强架构。
谢谢大家！！😁
使用 TensorFlow.js 进行时间序列预测
从在线 API 获取股票价格，并使用递归神经网络和长短期记忆(LSTM)以及 TensorFlow.js 框架进行预测
如今，机器学习变得越来越受欢迎，世界上越来越多的人将它视为一个神奇的水晶球:预测未来何时会发生什么。该实验使用人工神经网络来揭示股票市场趋势，并展示了时间序列预测根据过去的历史数据预测未来股票价格的能力。
免责声明:由于股票市场波动是动态的和不可预测的，由于多种因素，这个实验是 100%的教育，绝不是一个交易预测工具。
探索演示。
项目演练
本项目演练分为 4 个部分:
- 从在线 API 获取股票数据
- 计算给定时间窗口的简单移动平均值
- 训练 LSTM 神经网络
- 预测并将预测值与实际值进行比较
获取股票数据
在我们可以训练神经网络并做出任何预测之前，我们首先需要数据。我们正在寻找的数据类型是时间序列:按时间顺序排列的数字序列。获取这些数据的好地方是 alphavantage.co 的。这个 API 允许我们检索过去 20 年中特定公司股票价格的时序数据。
该 API 生成以下字段:
- 公开价格
- 那天的最高价
- 当天的最低价
- 收盘价(本项目中使用)
- 卷
为了给我们的神经网络准备训练数据集，我们将使用收盘股票价格。这也意味着我们将致力于预测未来的收盘价。下图显示了微软公司 20 年来的每周收盘价。
20 years of Microsoft Corporation weekly closing prices data from alphavantage.co
简单移动平均线
对于这个实验，我们使用监督学习，这意味着将数据输入到神经网络，它通过将输入数据映射到输出标签来学习。准备训练数据集的一种方法是从时序数据中提取移动平均值。
简单移动平均线(SMA) 是一种通过查看某个时间段内所有值的平均值来确定该时间段趋势方向的方法。时间窗口中的价格数量是通过实验选择的。
例如，假设过去 5 天的收盘价是 13，15，14，16，17，那么 SMA 就是(13+15+14+16+17)/5 = 15。因此，我们的训练数据集的输入是单个时间窗口内的一组价格，其标签是这些价格的计算移动平均值。
让我们计算微软公司每周收盘价数据的 SMA，窗口大小为 50。
这是我们得到的，蓝色的是每周股票收盘价，橙色的是 SMA。因为 SMA 是 50 周的均线，比周线价格平滑，周线价格可以波动。
Simple Moving Average of Microsoft Corporation closing prices data
培训用数据
我们可以用每周股票价格和计算出的 SMA 来准备训练数据。假设窗口大小为 50，这意味着我们将使用每连续 50 周的收盘价作为我们的训练特征(X)，并将这 50 周的 SMA 作为我们的训练标签(Y)。看起来是这样的…
接下来，我们将数据分成两组，训练集和验证集。如果 70%的数据用于训练，那么 30%用于验证。API 返回给我们大约 1000 周的数据，因此 700 周用于训练，300 周用于验证。
训练神经网络
既然训练数据已经准备好，是时候为时间序列预测创建一个模型了，为此我们将使用 TensorFlow.js 框架。TensorFlow.js 是一个用 JavaScript 开发和训练机器学习模型的库，我们可以在 web 浏览器中部署这些机器学习能力。
选择顺序模型，它简单地连接每一层，并在训练过程中将数据从输入传递到输出。为了让模型学习时序数据，创建了递归神经网络 (RNN)层，并在 RNN 中增加了多个 LSTM 细胞。
该模型将使用亚当 ( 研究论文)进行训练，这是一种流行的机器学习优化算法。均方根误差将决定预测值和实际值之间的差异，因此模型能够在训练过程中通过最小化误差进行学习。
下面是上面描述的模型的代码片段，Github 上的完整代码。
- 训练数据集大小(%):用于训练的数据量，剩余数据将用于验证
- 时期:数据集用于训练模型的次数(了解更多信息
- 学习率:每一步训练中的重量变化量(了解更多)
- 隐藏 LSTM 层:增加模型的复杂性，以便在更高维度的空间中学习(了解更多)
Web frontend, showing parameters available for tweaking
单击“开始培训模型”按钮…
User interface showing training model progress
该模型似乎在大约 15 个历元处收敛。
确认
既然模型已经定型，就该用它来预测未来值了，对于我们的情况，它是移动平均值。我们将使用 TFJS 中的 model.predict 函数。
数据被分成两组，训练组和验证组。训练集已用于训练模型，因此将使用验证集来验证模型。由于模型尚未看到验证数据集，如果模型能够预测接近真实值的值，那就太好了。
因此，让我们使用剩余的数据进行预测，这可以让我们看到我们的预测值与实际值有多接近。
The green line denotes the prediction of the validation data, from web demo
看起来模型预测的价格(绿线)与实际价格(蓝线)非常接近。这意味着该模型能够预测该模型看不到的最后 30%的数据。
可以应用其他算法，并使用均方根误差来比较 2 个或更多模型的性能。
预言；预测；预告
最后，模型已经过验证，预测值与其真实值非常接近，我们将使用它来预测未来。我们将应用相同的 model.predict 函数，并使用最后 50 个数据点作为输入，因为我们的窗口大小是 50。由于我们的训练数据每天都在增加，我们将使用过去 50 天作为输入，来预测第 51 天。
Predict the 51st day
为什么我的模特不表演了？
该模型过去从未出现过类似的数据。2020 年 3 月，市场下跌，并在一两个月内恢复，这在历史上从未发生过。该模型很可能无法预测这些时期股票价格的剧烈变化。
我们可以添加更多的功能。一般来说，更多的功能往往会使模型性能更好。我们可以包括交易指标，如移动平均线收敛背离(MACD)，相对强弱指数(RSI)，或布林线。
增加更多功能。Alpha Vantage API 提供的另一个惊人的 API 是基础数据。这意味着您还可以包括感兴趣的公司的年度和季度损益表和现金流。谁知道呢，这些功能可能会有用。
模型无法学习和预测的原因可能还有很多。这是机器学习的挑战；建立好的表演模型是一门艺术也是一门科学。
结论
除了使用简单的移动平均线，还有许多方法可以进行时间序列预测。未来可能的工作是利用来自各种来源的更多数据来实现这一点。
有了 TensorFlow.js，在网络浏览器上进行机器学习就成为可能，而且这实际上非常酷。
探索 Github 上的演示，这个实验是 100%教育性的，绝不是交易预测工具。
40%的吸尘器，40%的看门人，20%的算命师。
towardsdatascience.com](/data-scientist-the-dirtiest-job-of-the-21st-century-7f0c8215e845) 
Python 中使用动态时间弯曲的时间序列层次聚类
让我们考虑下面的任务:我们有一堆均匀分布的不同长度的时间序列。目标是通过定义数据中呈现的一般模式来对时间序列进行聚类。
在这里，我想介绍一种解决这个问题的方法。我们将使用层次聚类和 DTW 算法作为时间序列的比较度量。该解决方案在人力资源数据(员工历史得分)上运行良好。对于其他类型的时间序列，DTW 函数可能比其他指标更差，如 CID(复杂性不变距离)，MAE 或相关性。
我将跳过对层次聚类和 DTW 算法的理论解释，而专注于我为什么选择这样的组合:
- 层次聚类简单、灵活、可调(链接标准),并允许我们不聚类所有轨迹
- DTW 方法让我们可以比较不同长度的时间序列，根据我的经验，它非常适用于不频繁的时间序列
好了，我们开始吧！我们的进口:
import random
from copy import deepcopy
from scipy import interpolate
import numpy as np
from dtaidistance import dtwimport matplotlib.pyplot as pltfrom _plotly_future_ import v4_subplots
import plotly.graph_objects as go
from plotly.subplots import make_subplots
时间序列生成的一些参数和我们的阈值:
- 轨迹数量我们必须聚类的轨迹数量
- 最小轨迹长度和最大轨迹长度任何轨迹的长度下限和上限
- 门槛我们 DTW 的门槛
NUM_OF_TRAJECTORIES = 200
MIN_LEN_OF_TRAJECTORY = 16
MAX_LEN_OF_TRAJECTORY = 40THRESHOLD = 0.50
为简单起见，我们所有的轨迹将位于-1 和 1 之间。此外，我添加了一些平滑。
for item in range(NUM_OF_TRAJECTORIES):
length = random.choice(list(range(MIN_LEN_OF_TRAJECTORY, MAX_LEN_OF_TRAJECTORY + 1)))
tempTrajectory = np.random.randint(low=-100, high=100, size=int(length / 4)).astype(float) / 100
oldScale = np.arange(0, int(length / 4))
interpolationFunction = interpolate.interp1d(oldScale, tempTrajectory)
newScale = np.linspace(0, int(length / 4) - 1, length)
tempTrajectory = interpolationFunction(newScale)
trajectoriesSet[(str(item),)] = [tempTrajectory]
注意，所有轨迹都存储为列表类型的字典值(为了方便，我们将开始把它们组合成组)。出于同样的原因，轨迹的名称存储为元组。
我们的算法如下:
- 我们找到一对最接近的实体(轨迹-轨迹或轨迹-簇或簇-轨迹或簇-簇)
- 如果它们的距离低于阈值，则将它们分组到单个簇中
- 重复步骤 1
- 如果我们在步骤 2 失败或者我们得到一个大的集群，我们就停止我们的算法(所以我们所有的轨迹都进入其中——这意味着我们的阈值非常大)
算法的第一部分:
trajectories = deepcopy(trajectoriesSet)
distanceMatrixDictionary = {}iteration = 1
while True:
distanceMatrix = np.empty((len(trajectories), len(trajectories),))
distanceMatrix[:] = np.nan
for index1, (filter1, trajectory1) in enumerate(trajectories.items()):
tempArray = []
for index2, (filter2, trajectory2) in enumerate(trajectories.items()):
if index1 > index2:
continue
elif index1 == index2:
continue
else:
unionFilter = filter1 + filter2
sorted(unionFilter)
if unionFilter in distanceMatrixDictionary.keys():
distanceMatrix[index1][index2] = distanceMatrixDictionary.get(unionFilter)
continue
metric = []
for subItem1 in trajectory1:
for subItem2 in trajectory2:
metric.append(dtw.distance(subItem1, subItem2, psi=1))
metric = max(metric)
distanceMatrix[index1][index2] = metric
distanceMatrixDictionary[unionFilter] = metric
字典距离矩阵字典帮助我们保持已经计算好的距离。
Numpy 数组 distanceMatrix 在每一步开始时用 np.nan 填充。只需要保持索引对和计算的距离之间的表示。向 distanceMatrixDictionary 添加相同功能后，可能会被删除。
这部分代码允许我们比较所有可能的选项——轨迹-轨迹或轨迹-簇或簇-轨迹或簇-簇:
metric = []
for subItem1 in trajectory1:
for subItem2 in trajectory2:
metric.append(dtw.distance(subItem1, subItem2))metric = max(metric)
上面最后一行— metric = max(metric) —是称为“完全关联”的关联标准。它对我来说更好，但你可以尝试其他标准，甚至定制它。
好了，距离已经计算好了，让我们继续分组。
我们找到最低的距离和一对提供这个距离的指数。
这里为了简单起见，我们将只使用一对(第一对)。甚至，如果我们对于相同的距离有两个、三个或更多对，其余的将在下一次迭代中逐步处理。
minValue = np.min(list(distanceMatrixDictionary.values()))if minValue > THRESHOLD:
breakminIndices = np.where(distanceMatrix == minValue)
minIndices = list(zip(minIndices[0], minIndices[1]))minIndex = minIndices[0]
获得一对索引后，我们只需定义实体名称和值，组合它们，将组合放入字典中，并从字典中删除这些单个实体:
filter1 = list(trajectories.keys())[minIndex[0]]
filter2 = list(trajectories.keys())[minIndex[1]]trajectory1 = trajectories.get(filter1)
trajectory2 = trajectories.get(filter2)unionFilter = filter1 + filter2
sorted(unionFilter)trajectoryGroup = trajectory1 + trajectory2trajectories = {key: value for key, value in trajectories.items()
if all(value not in unionFilter for value in key)}distanceMatrixDictionary = {key: value for key, value in distanceMatrixDictionary.items()
if all(value not in unionFilter for value in key)}trajectories[unionFilter] = trajectoryGroup
之后，我们重复上一步，直到我们没有任何集群。
我已经描述了一般的方法，但是这个算法可以被简化、增强和修改以避免任何重新计算。
结果，我们得到这样的分组:
在这个集群中，我们看到 3 个不同长度的时间序列。它们都有相同的一般模式:前三分之一是局部最小值，后半部分是全局峰值，最后是全局最小值。
更多结果(此处，对于每个聚类，左边的子图表示原始轨迹长度，右边的子图——重新缩放到 MAX_LEN_OF_TRAJECTORY 用于比较):
Cluster #1–2 items
Cluster #2–2 items
Cluster #3–3 items
Cluster #4–3 items
Cluster #5–3 items
根据阈值的值，我们可以使我们的集群更大(更一般化)或更小(更详细)。
如果当前方法在另一个数据集上表现不佳，我们可以改进什么？
- 我们可以尝试用另一种距离度量来代替 DWT
- 我们可以尝试对时间序列进行额外的处理:缩放/重新缩放、平滑或去除异常值
- 我们可以尝试使用不同的阈值
- 我们可以尝试改变联系标准
找到最佳超参数后，可以重构代码并加快计算速度。
这里有可以使用的代码:
https://github . com/avchauzov/_ articles/blob/master/1.1 . trajectoriesclustering . ipynb
Python 中的时间序列—指数平滑和 ARIMA 过程
TL；DR:在这篇文章中，你将学习执行时间序列分析的基本步骤，以及趋势、平稳性、移动平均线等概念。您还将探索指数平滑方法，并学习如何在非平稳数据上拟合 ARIMA 模型。
时间序列随处可见
场景 1:你负责一家披萨配送中心，你想知道你的销售是否遵循特定的模式，因为你觉得每周六晚上你的订单数量都会增加……
场景 2:你的公司正在销售一种产品，而你负责预测，或者说预测，在未来某个特定时刻，这种产品所需的供应……
情况 3:您正在监控一个数据中心，并且想要检测任何异常情况，比如可能导致服务器停机的异常 CPU 使用情况。您跟踪 CPU 使用率的曲线，并想知道何时出现异常…
在每一种情况下，你都在处理时间序列。分析序列是一项迷人的工作，因为尽管有所有的数学模型(包括神经网络)，我们人类仍然无法预测未来，并且必须处理不确定性。让我们仔细看看什么是时间序列，可以用哪些方法来分析它们。在本文中，我们将广泛依赖于用 Python 编写的【stats models】库。
时间序列是按时间排序(或索引)的数据序列。它是离散的，每个点之间的间隔是常数。
系列的性质和类型
趋势 :数据的长期增减。这可以被看作是粗略穿过数据的斜率(不必是线性的)。
季节性 :当一个时间序列受到季节性因素(一天中的某个小时、一周、一月、一年等)的影响时，称之为季节性的。).季节性可以通过固定频率的良好循环模式来观察。
周期性 :当数据呈现非固定频率的上升和下降时，出现一个周期。这些波动通常是由经济条件造成的，而且往往与“商业周期”有关。这些波动的持续时间通常至少为两年。
残差 :每个时间序列可以分解成两部分:
-一个预测，由一个或几个预测值*
-残差。它们是每个时间步的观测值和预测值之间的差值。记住这一点*
t 时刻序列值= t 时刻预测值+t 时刻残差
**
Left: Series with a trend. Right: Series exhibiting seasonal behavior
**
Left: Cyclical series. Right: Random walk
时间序列的分解
每个时间序列都可以看作是几个部分的混合:
- 趋势(曲线长期向上或向下移动)
- 季节性成分
- 残差
下面是我们的时间序列分解成上述部分后的样子:
我们的时间序列显示出明显的上升趋势。它不是季节性的，你可以看到季节性的部分看起来很难看。残差的方差似乎随着时间的推移而增加，表明序列在结束时表现出更多的随机行为。
黄金法则:平稳
在我们进一步分析之前，我们的系列必须使 静止。
平稳性 是展现恒定统计性质(均值、方差、自相关等)的性质。).如果时间序列的平均值随着时间的推移而增加，那么它就不是稳定的。
用于静态化数据的转换:
- 去趋势 :我们去除序列中的潜在趋势*。这可以通过多种方式实现，具体取决于数据的性质:
- 指数化数据:以货币计量的数据与价格指数挂钩或与通货膨胀相关。因此，将序列除以该指数(即缩减)是去除数据趋势的解决方案。
- 非指数型数据:是否需要估计趋势是常数、线性还是指数型。前两种情况很容易，对于最后一种情况，有必要估计增长率(通货膨胀或通货紧缩),并应用与指数化数据相同的方法。*
- 差分 :可以通过减去周期值来去除季节性或周期性模式。如果数据是 12 个月的季节性数据，用 12 个滞后的差值序列减去该序列将得到一个“更平”的序列
- 日志 :在趋势中的复合率不是由于某个价格指数(即序列不是以货币计量的)的情况下，日志可以帮助线性化一个具有指数趋势的序列(回想一下 log(exp(x)) = x)。与通缩不同，它并没有消除任何最终趋势。
检查平稳性
绘制滚动统计数据
绘制滚动平均值和方差是直观检查我们系列的第一个好方法。如果滚动统计显示出明显的趋势(向上或向下)并显示出变化的方差(增加或减少幅度)，那么您可能会得出结论，该序列很可能不是平稳的。
扩充迪基-富勒试验
该检验用于评估时间序列是否是平稳的。不必深入假设检验的太多细节，你应该知道这个检验将给出一个叫做“检验统计”的结果，根据这个结果，你可以用不同的置信水平(或百分比)来判断时间序列是否平稳。
KPSS
KPSS(科维亚特科夫斯基-菲利普斯-施密特-申)检验测试的是序列是趋势平稳的零假设。换句话说，如果检验统计量的 p 值低于 X%置信度阈值，这意味着我们可以拒绝这个假设，并且序列不是具有 X%置信度的趋势平稳序列。高于阈值的 p 值将使我们接受这一假设，并得出序列是趋势平稳的结论。
自相关图(ACF 和 PACF)
自相关 (ACF)图表示具有自身滞后的序列的自相关。
偏自相关 (PACF)图表示一个序列和其自身滞后之间的相关量，该相关量不能用所有更低阶 - 滞后的相关来解释。理想情况下，我们不希望序列与其自身的滞后之间存在相关性。从图形上来说，我们希望所有的尖峰都落在蓝色区域。
正如我们所见，蓝色区域上方有几个尖峰，这意味着在滞后 1、2、3 和 4 处存在相关性。
选择模型
指数平滑法适用于非平稳数据(即具有趋势的数据和季节性数据)。
ARIMA 模型应该只用于静态数据。因此，应该去除数据的趋势(通过缩减或记录)，然后查看有差异的序列。
平滑方法
平滑方法作为加权平均值工作。预测是过去观察的加权平均值。权重可以是统一的(这是一个移动平均值)，或者遵循指数衰减-这意味着给予最近的观察更多的权重，给予旧的观察更少的权重。更高级的方法包括预测中的其他部分，如季节性成分和趋势成分。
我们将使用数学方程的分量形式。 y 表示我们的时间序列， p 表示我们的预测， l 表示水平， s 表示季节性成分， b 表示趋势成分
简单指数平滑
>什么时候使用？
数据点少，数据无规律，无季节性或趋势性。
*>数学背后 只要记住 SES 只有一个组件叫做 level (下面用“alpha”表示平滑参数)*。是前一级和当前观测值的加权平均值:
霍尔特线性平滑
>什么时候使用？ 数据中的趋势，无季节性。
数学背后 预测由一个水平分量和一个趋势分量组成。：
**
霍尔特阻尼趋势
霍尔特的线性趋势法的问题是，趋势在未来是恒定的，无限增加或减少。对于长期预测来说，这可能是个问题。因此，阻尼趋势法是一种添加阻尼参数的方法，以便趋势在未来收敛到恒定值(它使趋势变平)。参数𝑏由𝜙𝑏代替
>什么时候使用？
数据有趋势。使用乘法版本，除非之前已经记录了数据。在这种情况下，使用附加版本
Notice how the green curve flattens at the end of the forecast
霍尔特-温特季节平滑法
这将在下一篇文章的中讨论。但你要知道，这种方法就像霍尔特的线性平滑，除了我们增加了一个季节的成分！
ARIMA
ARIMA 模型(包括 ARMA、AR 和 MA 模型)是预测平稳时间序列的一类通用模型。ARIMA 模型由三部分组成:
- 序列((AR)部分)的滞后 值 的加权和**
- 序列(()MA***)的滞后 预测误差 的加权和***
- 一个 差 的时间序列((I)部分)
一个 ARIMA 模型通常被标注为【p，d，q】其中 p 表示 AR 部分、 d 差分的顺序(I部分)，而 q 表示 MA 项的顺序。
1)选择差分阶数
拟合 ARIMA 模型的第一步是确定差分阶数以平稳化序列。为此，我们来看看 ACF 和 PACF 图，并记住这两条规则:
”—规则 1:如果序列对大量滞后具有正自相关，那么它可能需要更高阶的差分。
—规则 2:如果 lag-1 自相关为零或负值，或者自相关都很小且无模式，则序列不需要高阶差分。如果滞后-1 自相关为-0.5 或更负，则该系列可能有过度差异。 (罗伯特·诺，统计预测)
我们从记录数据开始，因为原始数据呈现指数趋势:
记录的序列似乎更平坦，但它是平稳的吗？让我们计算一个 KPSS 测试来检验这一点:
检验统计量高于临界值，我们拒绝零假设，我们的序列是而不是趋势平稳的。然而，为了使这篇文章简短，我们将继续好像它是。
(—为了得到一个趋势稳定的序列，我们可以考虑线性回归我们的对数序列，并用回归系数除我们的序列…— )
现在让我们来看看记录的一阶差分数据的 ACF 图:
“ACF —记录的数据”图表显示了非平稳数据，其特征是峰值缓慢线性衰减(参见上述规则 1)。添加一阶差会在滞后值 1 处产生单个负尖峰。根据规则 2，我们不需要进一步区分这个系列。让我们通过比较(0，0，0)和(0，1，0) ARIMA 模型来检查我们的结果:
****
ARIMA(0，1，0)的阿凯克信息标准(AIC)较低，这意味着该模型的表现优于 ARIMA(0，0，0)。让我们看看残差并检查它们的方差:
Which residuals do you think look better ?
2)选择移动授权订单
现在我们知道我们必须在模型中包含一阶差分，我们需要选择移动平均阶。这是通过查看差分序列来完成的(因为我们刚刚看到一阶差分序列是平稳的)。同样，我们看看我们的 ACF 和 PACF 图，记住这条规则:
“如果差分序列 ACF 的滞后-1 自相关为负，和/或存在锐截止，则选择顺序为 1 的 MA ”。
问:为什么选择 1 阶而不是更高阶的移动授权？因为我们会一步一步来。如果我们在更高的滞后下观察自相关，并且通过查看我们的(1，1，0)模型的自相关残差，我们仍然观察到这些尖峰，我们可以增加我们的 MA 阶数，尽管通常不建议超过 2！
注意 AIC 是如何再次下降的，残差方差是如何减小的。这表明我们的(1，1，0) ARIMA 比(0，1，0)模型表现得更好！
3)选择应收账款订单
现在你可能会想:我们有必要添加一个 AR 术语吗？答案是否定的。事实上，在以下情况下，您应该添加一个 AR 术语:
如果差分序列 PACF 的滞后-1 自相关为负，和/或有一个锐截止，那么选择一个为 1 的 AR 阶。
在我们的例子中，我们在 ACF 和 PACF 中都观察到负的滞后-1 自相关。请注意，差分序列的 PACF 在滞后 1 和滞后 2 处显示了两个负峰值，这意味着理论上我们可以将 AR 阶数提高到 2。
这是我们通过拟合(1，1，1) ARIMA 得到的结果:
****
添加 AR 项降低了 AIC，方差从 0.155 降到了 0.145，我们做得很好！我们应该增加一个 AR 项并选择一个(2，1，1) ARIMA 吗？让我们来看看 ARIMA(1，1，1)残差的 ACF & PACF:
图中没有显著的自相关滞后值，我们的模型似乎不需要任何额外的 AR 项。你可能会指出第 12 阶段的峰值:这可能是季节性的迹象。我们将在的下一部分中与萨里玛和萨里马克斯一起探讨这个问题。
绘制预测
幸运的是，statsmodels 有一个很好的 API，允许我们从 ARIMA 过程中绘制预测。我强烈建议您将数据放在带有 DateTimeIndex 的 DataFrame 中，因为 plot_predict() 方法确实喜欢日期…
在这种情况下，我们将选择 12 个预测步骤，并将 dynamic 关键字设置为 True:这将强制算法使用其自己的预测值作为未来预测的滞后值:
让我们看看更大的预测窗口(34 个预测步骤):
那都是乡亲们！在下一部分的中，我将讨论季节性 Holt-Winter 和季节性 ARIMA。希望你喜欢:)
附加讲座:
如果你对时间序列预测和分析感兴趣，我强烈推荐下面两个:
- 预测、原理与实践(Hyndman&Athanasopoulos)
- 【罗伯特·诺】****
你也可以查:
如果你想知道 2019 年哪支球队真的是最好的 NBA 球队，请查看 这篇文章 ！
Python 中的时间序列—第 2 部分:处理季节性数据
在 第一部分 中，你学习了趋势和季节性、平滑模型和 ARIMA 过程。在这一部分中，您将学习如何处理季节性模型，以及如何实现季节性 Holt-Winters 和季节性 ARIMA (SARIMA)。
获取数据
我们将使用“每月牛奶产量”数据:
季节分解
在上一篇中，我简单讲过季节分解。季节性分解的思想是，任何序列都可以分解为三个部分的总和(或乘积):趋势、季节性成分和残差。
在我们的例子中，我们将使用 statsmodels 提供的季节性分解函数:
趋势分量是一条增长曲线，似乎达到了一个平稳状态，并最终在终点下降。这是一个有用的信息，当我们将消除趋势，以稳定的数据。
将数据固定化
应对趋势
基于我们的分解，我们看到趋势是跟随一个向上的运动，然后在最后达到稳定状态(或者下降)。为了消除这种趋势，我们将尝试一种新颖的方法，即回归 STL 分解给出的趋势。然后，如果回归令人满意，我们将尝试通过从原始序列中减去获得的回归来“缩小”序列。
我们回归的 r 平方系数相当好(0.988)，但是查看残差，我们注意到它们的方差从左到右增加。这是异方差的标志。您可以自己查找确切的定义，只要记住这并不好，我们的模型中缺少了一些东西。我们将放弃回归趋势的想法，因为它现在需要太多的努力。
应对趋势，第 2 页
我们的系列仍然需要站化，我们将回到基本方法，看看我们是否可以消除这种趋势。我们将分几个步骤操作:
- 对原始数据进行 ADF 测试以检查平稳性
- 12 个月差异的 ADF 测试
- 对记录数据的 12 个月差异进行 ADF 检验
- 对减去 12 个月移动平均值的数据进行 ADF 测试
如果您需要代码，请参考第一部分，了解如何进行 ADF 测试
让我们用 KPSS 检验来确认数据减去 12 个月的均线是平稳的:
测试统计低于阈值，因此我们的序列是平稳的(回想一下，对于 KPSS 测试，零假设是序列是平稳的)。
霍尔特-温特季节平滑模型的拟合
记得从第一部分得知，霍尔特-温特的模型有几个部分:一个 水平 ，一个 趋势 ，在季节性平滑的情况下，一个 季节性 成分。这背后的数学有点难，所以我不会把它放在这里，只要记住上面的三个组成部分。以下是我们拟合模型后得到的结果:
没有阻尼的模型似乎表现更好，因为 AIC 更低。
季节性 ARIMA
老 ARIMA 又回来了，但这一次我们将为它添加季节性元素！SARIMA 模型由两部分组成:非季节性部分和季节性部分。如果你还记得上一部分，非季节性部分由 AR§项、MA(q)项、I(d)项组成，我们写 ARIMA (p，d，q)。季节性的部分也是一样的！
我们的舍利玛是指 ARIMA (p，D，q)x(P，D，Q)。
问题:我们如何确定所有这些 P，D，Q，P，D，Q？
每当你听到“ARIMA”，就想到“ACF/PACF”。
有了 ACF 和 PACF 图，我们将能够猜测我们的参数的合理值。让我们绘制我们的静态化数据的 ACF 和 PACF 图:
鉴于在 ACF 中观察到的尖峰信号的缓慢衰减，数据显然不是稳定的。这是一个信号，表明我们应该对 x: 取一个额外的差值
我们的系列现在是静止的。总而言之，我们有:
- 取平稳序列的 12 个滞后差。这相当于“一个季节滞后”的差异。
- 取 12 滞后差分序列的 1 滞后差分。这对应于 1 的非季节性差序。
这给了我们一个萨里玛(？, 1, ?)x(？, 1, ?).
AR 和 MA 术语呢？
嗯，我们知道在 ACF 和 PACF 中仍然有显著的峰值，所以我们可以决定添加 AR 术语和 MA 术语。我们应该从哪一个开始？这一次，我们将使用网格搜索 AR 和 MA 术语的可能组合(最大顺序限制为 2)…
当运行不同的配置时，使用 SARIMA(0，1，2)x(0，1，2)可以获得最低的 AIC。
使用*plot _ diagnotis()*命令，我们可以显示残差并检查它们的正态性(记住，残差应该是正态分布的(零均值和一元方差)并且不相关)。
请注意正态 Q-Q 图和直方图如何显示我们的残差遵循正态分布(我们可以使用额外的测试来测试正态性，如 Ljung-Box)。相关图还显示残差是不相关的！
附加讲座:
如果你对时间序列预测和分析感兴趣，我强烈推荐下面两个:
- 预测、原理与实践 (海德曼& Athanasopoulos)
- (罗伯特·诺)
你还可以查看 这篇文章 关于神经图灵机和 这篇文章 关于基于注意力的神经网络！
Python 中的时间序列第 3 部分:用 LSTMs 预测出租车出行
在本文中，我们将使用 LSTMs 预测纽约的 FHV 数量(如优步出租汽车)。我们将学习如何使用蒙特卡洛法将置信区间添加到我们的预测中。
Photo by Lexi Anderson on Unsplash
介绍
LSTM(Long Short Memory)是一种 a 型递归神经网络(RNN)架构，由 Sepp Hochreiter 和 Jürgen Schmidhuber 于 1997 年提出。rnn 是专门设计用于通过递归机制处理顺序数据的深度神经网络。它们以一种自回归的方式表现，因为它们通过内部状态来跟踪过去(因此是“记忆”部分)。它们被广泛用于语音识别、机器翻译、语音合成等。但是当用于时间序列时，LSTMs 值多少呢？假设可用数据的规模足够大，它们可以被证明对非线性关系的建模非常有用…
优步用例:贝叶斯预测
在找用 LSTMs 实现时间序列预测的论文时，找到了优步 2017 年写的一篇论文，优步对时间序列的深度自信预测。这篇论文背后的基本问题是:我们有多大信心(即我们如何量化不确定性)用 LSTMs 进行预测？
优步开发的方法是编码器-解码器(用作自动编码器)和完全连接的前馈网络的混合，用于根据以前的数据预测城市的出行次数，或实时检测异常。
The Uber LSTM forecasting architecture (Zhu & Laptev, 2017)
优步的论文是最先使用贝叶斯方法进行时间序列预测的论文之一。如果你想了解更多关于贝叶斯神经网络和贝叶斯推理的知识，可以看看下面的链接:
- 让你的神经网络说不知道
- 辍学为贝叶斯近似
- 深度贝叶斯神经网络
- 黑客的贝叶斯方法，卡梅隆·戴维森-皮隆
贝叶斯神经网络
简而言之，贝叶斯神经网络估计每个权重的概率分布，而经典神经网络试图找到每个权重的最佳值。听到“贝叶斯”，就想到“概率分布”。
——但是预测不确定性和贝叶斯网络有什么联系呢？
想象两个气象专家，鲍勃和爱丽丝。你们都知道他们很擅长预测天气温度，但有时他们也会弄错。你问他们明天早上 8 点的温度是多少，这样你就可以知道你是否需要穿上外套。
鲍勃说:“将会是 15.6 摄氏度”。爱丽丝说:“我有 95%的把握气温会在 16 到 18 摄氏度之间”。虽然他们似乎不同意，但你会相信谁呢？
就我个人而言，我相信爱丽丝有两个原因:
- 她给了我一个信心区间。当有人能够告诉我一些事情以及他/她对这些信息有多自信时，我会感到更加放心。
- 我并不在乎确切的温度，因为 1 摄氏度的温差不会影响我穿上外套的决定。
然而，如果 Alice 告诉我“我 95%确定温度将在 0°C 到 18°C 之间，我会说虽然她给了我一个信心水平，但这个区间太大了，无法提供信息……
BNN 框架下的不确定性
我们通常将不确定性分为 3 类: 模型不确定性 ， 固有噪声 ， 模型误设定 。前两个是最著名的，通常被称为和 任意 不确定性。
模型(或认知*)不确定性是关于我们的模型参数的不确定性。数据越多，就越能解释(即“减少”)这种不确定性。噪声(或任意性)不确定性是指观测值中的噪声。如果所有样本的噪声不确定度都是常数，我们称之为 同方差 不确定度。否则，如果一些样本比其他样本更不确定，我们将使用术语 异方差随机变量 。噪声不确定性不能随着更多的数据而减少。*
在优步的论文中，使用了第三种不确定性:模型设定错误。这种不确定性旨在“在预测与训练数据集有很大不同的模式的未知样本时捕捉不确定性”。
— 现在为什么要区分所有这些不确定性 ？
嗯，正是因为有些可以用更多的数据(模型/认知)来对抗，而有些不能。因此，一些研究人员认为，考虑到即使有更多的数据也无法减少随机不确定性，关注随机不确定性更有意义( Kendall & Gal，2017 )
在本文中，为了简单起见，我们将只关注模型的不确定性。优步使用不同的算法来评估另外两种不确定性，但是研究它们超出了本文的范围。
获取和准备数据
论文中的作者使用美国 8 个城市 4 年的数据来训练他们的模型。3 年用于培训，1 年用于测试。不幸的是，优步还没有发布这个数据，但是为了重现他们论文的结果，我们将使用纽约开放数据门户上的可用数据。我们选取了三年的数据，时间跨度从 2015 年初到 2018 年年中。每天对数据进行重新采样。这是完整数据的样子
平稳性
如果你习惯于分析时间序列，我们可以注意到一些你应该熟悉的东西:
- 我们观察到明显的上升趋势
- 均值和方差随时间增加
- 我们还观察到可能由外部事件(节假日和天气)引起的峰值。)
前两个要点表明我们的系列显然不是一成不变的。后者表明我们可能需要将外部数据合并到我们的系列中。平稳性可以用扩展的迪基-富勒测试或 KPSS 测试来检验。那么我们应该像对 ARIMA 模型那样进行这些测试吗？
—答案是:不一定。
当我们平稳化一个序列时，我们想要的是在整个时间内均值和方差不变。这保证了如果你取一个 N 点样本进行预测，并用另一个 N 点不同的样本重复这一过程，那么你的 N 点样本和你试图预测的 N+1 点之间的关系是相同的(即它们来自相同的分布)。如果均值和方差不相等，那么您是从不同的分布中抽取样本来进行预测，您的模型肯定无法进行概括。
与 ARIMA 不同，rnn 能够模拟数据中的非线性关系。
rnn，尤其是 LSTM 和 GRU，能够捕捉长期依赖关系(前提是你有足够多的数据！).因此，去趋势化和去季节性之类的问题不那么重要，但你应该经常问自己:
我用于测试的数据是否遵循与我用于训练的数据相同的行为(即来自相同的分布)？
假期和天气
添加假日指示非常简单。我们使用 Python 中的 holidays 库来获取 2015 年到 2018 年的假期日期。
在我们的系列中添加了一个假日布尔之后，我们仍然观察到无法解释的峰值…
在互联网上快速浏览显示，其中几天实际上是纽约州遭受暴风雪等极端天气事件袭击的日子。
我们在下面的图表中用红色标出了节假日，绿色标出了极端天气事件:
因此，在我们的数据框架中，我们有一个“计数”列，一个“是 _ 假日”列，一个“是 _ 恶劣天气”列。但是，假设我们要进行预测，我们需要“预测”这些日期，就像我们对未来预测一样，因此我们将创建两个额外的列，指示第二天是假日或第二天预计会出现极端天气:
*weather = [datetime.datetime.strptime(date, "%Y-%m-%d") for date in ['2018-01-04', '2018-03-21','2017-03-14','2017-02-09','2016-01-23']]holidays = [date for y in range(2015, 2019) for date, _ in sorted(holidays.US(years=y).items())]df['is_holiday'] = np.where(df.index.isin(holidays), 1, 0)
df['bad_weather'] = np.where(df.index.isin(weather), 1, 0)
df['next_is_holiday'] = df.is_holiday.shift(-1)
df['next_bad_weather'] = df.bad_weather.shift(-1)*
选择窗口大小和回测
做时间序列预测时，你可能会听说 回测 。回溯测试是在训练过程中使用的一个程序，它包括以递增的方式将数据分割成块。在每一次迭代中，一个块被用作你的训练集。然后你试着在你的块之前预测一个或多个值。可以使用两种方法，展开和滑动窗口:
source: https://www.datapred.com/blog/the-basics-of-backtesting
在我们的案例研究中，作者使用由步长等于 1 的 28 天滑动窗口组成的样本，用于预测下一个值(提前 1 步预测)。
记录和缩放
在论文中，作者首先采用数据的对数来“减轻指数效应”。然后，他们在每个窗口中减去窗口的第一个值，以消除趋势，并根据窗口第一个值的波动来训练网络(等式(1) )。也可以想到其他方法，例如减去第一值并除以第一值(等式(2) ):
构建数据集
我们的数据集将是一个生成器，产生一批滑动窗口(每一批都是在未来移动 1 个值的前一批)。为了遵循论文的说明，我们还将在每一批中从所有其他值中减去第一个值。然后，每批样本被分为 28 天样本和 1 天目标样本。
定义模型
我们将使用 PyTorch 来定义我们的模型。一个简单的原因是，我们将在推理过程中使用 dropout，并且它很容易在 PyTorch 中实现。我们将从使用一个简单的 LSTM 网络开始，如论文中所定义的:1 个具有 128 个单元的 LSTM 层，1 个具有 32 个单元的 LSTM 层，以及一个具有 1 个输出的全连接层。在每一层之后添加漏失。
——中途辍学者 中途辍学可以被看作是做贝叶斯推理的一种方式(尽管围绕这一点还有争论)。从技术上讲，退出是神经网络中使用的过程，包括随机退出单元(以及它们的连接)。
随机打开和关闭神经元的事实大致相当于执行伯努利分布的采样，因此“模拟”了贝叶斯神经网络的机制(其中权重是分布，而不是单个值)。应用 dropout 有点像我们从网络中“取样”。如果我们在推理过程中多次重复这个过程，我们将得到不同的预测，用这些预测我们可以估计一个分布，最终，不确定性！总而言之:
现在，让我们通过在每个 LSTM 层之间添加下降层来定义模型(注意 train 如何设置为 True，以便在训练和测试期间使用下降)
培养
我们用 Adam 优化器和学习率设置为 0.001，batch_size 为 1，训练 5 个时期。
结果:
Fitted values on the train set
测试— 1 天预测范围
本文的主要兴趣之一是不确定性的估计。
具体而言，在每个隐藏层之后应用随机漏失，并且模型输出可以近似地视为从后验预测分布生成的随机样本。因此，可以通过几次重复的模型预测的样本方差来估计模型的不确定性。
因此，本文背后的想法是用随机漏失运行几次模型，这将产生不同的输出值。然后，我们可以计算输出的经验均值和方差，以获得每个时间步长的置信区间！
对于每一步，我们预测 100 个值。所有的值都是不同的，假设我们保持辍学集。这使我们能够模拟从我们的网络中采样(事实上，辍学与贝叶斯神经网络密切相关，因为通过随机断开权重，它模拟了一种概率分布)。我们将这 100 个值的平均值作为我们的预测平均值，并将这 100 个值的标准偏差用于我们的置信区间。假设我们的预测来自正态分布 N ( μ ，σ)——均值 μ 等于我们的经验均值，标准差 σ 等于我们的经验标准差——那么我们可以估计置信区间。在我们的例子中，它由下式给出:
**
Predicted values and uncertainty intervals
预测和测试曲线似乎非常接近！然而，我们需要找到一个度量标准来查看我们的模型表现如何。
如果我们查看 95%预测区间的经验覆盖率(即预测的 95%置信区间中包含的真实测试值的数量)，我们会获得 28.51%的值，与论文中测试集上获得的值相差甚远……当我们取 99% CI 时，该值会上升到 44%覆盖率。这不是一个很好的结果，但请记住，我们的置信区间很小，因为我们只预测模型的不确定性…
测试— 7 天预测范围
单日预测的结果对于预测模型来说并不太好，而且我们对模型的要求并不高，因为我们要求的是一个非常短的预测范围。如果我们训练我们的模型来预测更大范围的预测会怎么样？
结论
lstm 因其捕捉长时间相关性的能力而显示出建模时间序列的有趣前景，但应仅用于大型数据集。永远不要期望你的 RNN 在 100 个样本的数据集上给出好的结果！神经网络中的参数数量和可用数据大小之间的平衡对于避免过度拟合非常重要。优步论文的作者之一尼古拉·拉普捷夫总结说:
经典模型最适合:
○短的或不相关的时间序列
○已知的世界状态神经网络最适合:
○大量时间序列
○长时间序列
○隐藏交互
○解释不重要
时间序列机器学习分析和需求预测与 H2O 和 TSstudio
我如何使用机器学习来实现每周收入的时间序列预测
Photo by Ben Elwood on Unsplash
传统的时间序列分析和预测方法，如、 霍尔特-温特斯指数平滑 、ARMA/ARIMA/萨里玛 和 ARCH/GARCH ，已经建立了几十年，并在商业和金融**(例如预测股票价格和分析金融市场的趋势)**
最近，开源框架的普及和更广泛的可用性，如 Keras 、 TensorFlow 和 scikit-learn 帮助机器学习方法，如 【随机森林】 、 极端梯度提升、 时间延迟神经网络 这些技术允许通过一组时间延迟将历史信息作为输入引入模型。
与更传统的方法相比，使用机器学习模型的优势在于，它们可以有higher predictive power，尤其是当预测器与响应有明确的因果联系时。此外，它们可以处理大量输入的复杂计算。
不过，它们往往有wider array of tuning parameters，一般都是more complex比【经典】型号多，而且可以expensive来上手，无论是计算能力还是时间。更重要的是，他们的black box本性使得他们的输出更难解释，并催生了不断增长的 机器学习可解释性 (我不打算触及这一点，因为这超出了项目的范围)
项目结构
在这个项目中，我将详细解释用机器学习模型对时间序列数据建模所需的各个步骤。
其中包括:
- 探索性时间序列分析
- 特征工程
- 模型训练和验证
- 模型性能对比及预测
特别是，我使用TSstudio来执行“传统的”时间序列探索性分析来描述时间序列及其组成部分，并展示如何使用我收集的洞察力来为机器学习管道创建特征以最终生成每周收入预测。
对于建模和预测，我选择了高性能、开源的机器学习库。我正在安装各种机器学习模型，例如广义线性模型、梯度推进机器和随机森林，并且还使用 AutoML 进行自动机器学习，这是H2O库最令人兴奋的特性之一。**
数据
**library(tidyverse)
library(lubridate)
library(readr)
library(TSstudio)
library(scales)
library(plotly)
library(h2o)
library(vip)
library(gridExtra)**
我在这里使用的数据集附带了一个红皮书出版物，可以在附加资料部分免费下载。这些数据涵盖了样本户外公司的销售额orders，这是一家虚构的 B2B 户外设备零售企业，并提供了他们销售的products以及他们的客户(在他们的案例中是retailers)的详细信息。由于它的人工性质，这个系列呈现了一些古怪和怪癖，我将在整个项目中指出这些。
在这里，我只是加载已编译的数据集，但如果你想继续，我也写了一篇名为加载、合并和连接数据集的文章，其中我展示了我如何组装各种数据馈送，并整理出变量命名、新功能创建和一些常规内务任务等内容。
**# Import orders
orders_tmp <-
read_rds("orders_tbl.rds")**
你可以在我的 Github 库上找到完整的代码。
初步探索
时间序列数据有一组独特的特征，如时间戳、频率和周期/时期，它们在描述性和预测性分析中都有应用。r 提供了几个类来表示时间序列对象(xts和zoo仅举几个主要的例子)，但是为了涵盖这个项目的描述性分析，我选择了使用ts类，它由TSstudio库支持。****
TSstudio附带了一些非常有用的功能，用于时间序列对象的交互式可视化。我真的很喜欢这个库使用作为它的可视化引擎！**
首先，我选择我的分析所需的数据(在本例中是order_date和revenue，并将其聚合为一个周频率。我用一个_tmp后缀来标记这个数据集，表示它是一个临时版本，在它可以使用之前还需要一些步骤。**
**revenue_tmp <- orders_tmp %>%
# filter out final month of the series, which is incomplete filter(order_date <= "2007-06-25") %>%
select(order_date, revenue) %>%
mutate(order_date =
floor_date(order_date,
unit = 'week',
# setting up week commencing Monday
week_start = getOption("lubridate.week.start", 1))
) %>%
group_by(order_date) %>%
summarise(revenue = sum(revenue)) %>%
ungroup()revenue_tmp %>% str() ## Classes 'tbl_df', 'tbl' and 'data.frame': 89 obs. of 2 variables:
## $ order_date: Date, format: "2004-01-12" "2004-01-19" ...
## $ revenue : num 58814754 13926869 55440318 17802526 52553592 ...**
这个系列跨越了3&1/2 年的销售额orders，所以我应该期望至少有 182 个数据点，但是只有 89 个观察值！
让我们仔细看看发生了什么:
**revenue_tmp %>% head(10) ## # A tibble: 10 x 2
## order_date revenue
## <date> <dbl>
## 1 2004-01-12 58814754\.
## 2 2004-01-19 13926869\.
## 3 2004-02-09 55440318\.
## 4 2004-02-16 17802526\.
## 5 2004-03-08 52553592\.
## 6 2004-03-15 23166647\.
## 7 2004-04-12 39550528\.
## 8 2004-04-19 25727831\.
## 9 2004-05-10 41272154\.
## 10 2004-05-17 33423065.**
****这里需要注意一些事情:这个系列呈现出一种不同寻常的每周模式，平均每月记录两次销售。此外，2004 年的第一周没有任何销售记录。
为了使用ts对象进行时间序列分析，我需要确保我每年有整整 52 周的时间，这也应该包括没有销售额的几周。
所以在将我的数据转换成一个ts对象之前，我需要:
- 在系列开始处增加 1 次观察以确保第一年包括 52 次每周观察
- ****一切都是按时间顺序排列的。这一点尤其重要，因为输出可能无法正确映射到序列的实际索引，从而导致不准确的结果。
- 填补不完整日期时间变量的空白。
padr库中的pad函数为每个缺失的时间点插入一条记录(默认填充值为 NA) - 用零替换缺失值，因为空周没有销售记录。来自
padr库的fill_by_value函数对此有所帮助。
**revenue_tbl <-
revenue_tmp %>%
rbind(list('2004-01-05', NA, NA)) %>%
arrange(order_date) %>%
padr::pad() %>%
padr::fill_by_value(value = 0)revenue_tbl %>% summary() ## order_date revenue
## Min. :2004-01-05 Min. : 0
## 1st Qu.:2004-11-15 1st Qu.: 0
## Median :2005-09-26 Median : 0
## Mean :2005-09-26 Mean :24974220
## 3rd Qu.:2006-08-07 3rd Qu.:52521082
## Max. :2007-06-18 Max. :93727081**
现在我可以看看周收入，我的反应变量
**revenue_tbl %>%
ggplot(aes(x = order_date, y = revenue)) +
geom_line(colour = 'darkblue') +
theme_light() +
scale_y_continuous(labels = scales::dollar_format(scale = 1e-6,
suffix = "m")) +
labs(title = 'Weekly Revenue - 2004 to June 2007',
x = "",
y = 'Revenue ($m)')**
如前所述，该系列是人为生成的，并不一定反映现实生活中会发生什么。如果这是一个实际的分析咨询项目，我肯定会质疑我的客户每周的销售频率。
但假设这是真实的交易，这里的挑战是构建一个经过深思熟虑选择的有意义的特征，并在几个机器学习模型上测试它们，以找到一个能够产生良好预测的模型。
挑战除外！
探索性分析
在这一部分，我将探索时间序列，检查其组成部分和季节性结构，并进行相关性分析，以确定该序列的主要特征。
在将我的系列转换成一个ts对象之前，我需要定义系列的start(或end)参数。我希望从一年的第一周开始计算周数，这样可以确保所有事情都符合ts框架。
**start_point_wk <- c(1,1) start_point_wk
## [1] 1 1**
我通过选择响应变量(revenue)作为数据参数并指定 52 周的频率来创建ts对象。
**ts_weekly <-
revenue_tbl %>%
select(revenue) %>%
ts(start = start_point_wk, frequency = 52)ts_info(ts_weekly) ## The ts_weekly series is a ts object with 1 variable and 181 observations
## Frequency: 52
## Start time: 1 1
## End time: 4 25**
用ts_info()函数检查系列属性显示，该系列是一个包含 1 个变量和 181 个观察值的每周ts对象。
时间序列组件
现在让我们借助TSstudio的图形功能来绘制我们的时间序列。
ts_decompose将序列分解为其元素:趋势、季节性和随机成分。
**ts_decompose(ts_weekly, type = 'additive')**
- 趋势:该系列似乎没有周期性成分，但显示出明显的上升趋势。趋势可能不是线性的，我将通过在特性中包含一个平方趋势元素来捕捉这一点。****
- 季节性:这个情节显示了一个明显的季节性模式，我将在接下来探究。**
- 随机:随机成分看起来是随机分布的。**
季节性成分
现在让我们放大这个系列的季节性部分
**ts_seasonal(ts_weekly, type = 'normal')**
尽管有 12 个不同的“峰值”(一年中的每个月一个)，但该图并不表明存在典型的季节性。然而，除了极少数例外，销售记录在每年的同一周，我将尝试用一个特征变量来捕捉这种规律性。
相关分析
自相关函数(ACF) 描述了序列与其滞后之间的相关程度。
由于手头的系列的奇怪性质，AC 图不是很容易阅读和解释:它表明存在滞后结构，但由于系列中的噪声，很难拾取有意义的模式。
**ts_acf(ts_weekly, lag.max = 52)**
然而，我可以利用滞后可视化和玩滞后数字，以确定序列和它的滞后之间的潜在相关性。
在这种情况下，将周数与季度频率对齐显示了与季度滞后的明显的线性关系。为简单起见，我将只在模型中包括滞后 13 ，以控制上一季度销售水平的影响。
**ts_lags(ts_weekly, lags = c(13, 26, 39, 52))**
使用相同的技巧揭示了与第一年滞后的强线性关系。同样，我将在模型中只包括一个滞后 52 。
ts_lags(ts_weekly, lags = c(52, 104, 156))
探索性分析总结
- 该系列的购买频率为两周一次，两周一次，没有典型的季节性。然而，除了极少数的例外情况，每年的销售额都大致记录在同一个星期。****
- 该系列似乎没有周期性成分，但显示出明显的上升趋势以及潜在的非线性趋势。
- 由于数据嘈杂，ACF 很难解释，但滞后图暗示了年度和季度滞后结构。****
系统模型化
****建模和预测策略是:
- 使用 Q1 2007 训练并交叉验证 Q1 2007 之前的所有车型，并比较其样本内预测性能**。**
- 假设我没有 Q2 2007 的数据，使用所有合适的模型生成该期间的预测，并将它们的预测性能与 Q2 2007 实际值进行比较
下面是该策略的可视化表示。我从经验中发现，用一个好的形象化来支持解释是一个让你的观点变得生动的好方法，尤其是在进行时间序列分析的时候。
所有模型精度将与性能指标和实际与预测图进行比较。
我将使用的绩效指标是:
- R 是一个拟合优度指标，它以百分比的形式解释了由于特征变量的变化而导致的响应变量的变化量。**
- RMSE (或均方根误差)是残差的标准差，衡量预测误差的平均大小。基本上，它告诉你残差是如何分布的。**
revenue_tbl %>%
filter(order_date >= "2005-01-03") %>%
ggplot(aes(order_date, revenue)) +
geom_line(colour = 'black', size = 0.7) +
geom_point(colour = 'black', size = 0.7) +
geom_smooth(se = FALSE,
colour = 'red',
size = 1,
linetype = "dashed") +
theme_light() +
scale_y_continuous(limits = c(0, 11.5e7),
labels = scales::dollar_format(scale = 1e-6,
suffix = "m")) +
labs(title = 'Weekly Revenue - 2005 to June 2007',
subtitle = 'Train, Test and Forecast Data Portions',
x = "", y = 'Revenue ($m)') + # Train Portion
annotate(x = ymd('2005-12-01'), y = (10.5e7), fill = 'black',
'text', label = 'Train\nPortion', size = 2.8) + # Test Portion
annotate(x = ymd('2007-02-05'), y = (10.5e7),
'text', label = 'Test\nPortion', size = 2.8) +
geom_rect(xmin = as.numeric(ymd('2006-12-18')),
xmax = as.numeric(ymd('2007-03-25')),
ymin = -Inf, ymax = Inf, alpha = 0.005,
fill = 'darkturquoise') + # Forecast Portion
annotate(x = ymd('2007-05-13'), y = (10.5e7),
'text', label = 'Forecast\nPortion', size = 2.8) +
geom_rect(xmin = as.numeric(ymd('2007-03-26')),
xmax = as.numeric(ymd('2007-07-01')),
ymin = -Inf, ymax = Inf, alpha = 0.01,
fill = 'cornflowerblue')
****非常重要的一点:如你所见，我没有在图中显示 2004 年的数据。这是因为，每当您在模型中包含滞后变量时，用于计算滞后的第一个周期【drop off】数据集，将无法用于建模。在有年度滞后的情况下，所有观察值都提前一年，并且由于没有 2003 年的销售记录，这导致前 52 周被从分析中删除。
特征创建
现在，我可以开始将时间序列探索的结果整合到我的特征变量中。为此，我创造了:
趋势**特征:(一个趋势和趋势平方)。这是通过一个简单的数字指数来控制上升趋势和潜在的非线性趋势。**
Lag 特性:(a _lag 13 和 _lag 52 )捕捉观察到的收入与其季度和年度季节性滞后的相关性。
季节性**特性，用于处理两周工作、两周休息的采购频率
model_data_tbl <-
revenue_tbl %>%
mutate(trend = 1:nrow(revenue_tbl),
trend_sqr = trend^2,
rev_lag_13 = lag(revenue, n = 13),
rev_lag_52 = lag(revenue, n = 52),
season = case_when(revenue == 0 ~ 0,
TRUE ~ 1)
) %>%
filter(!is.na(rev_lag_52))
下一步是创建训练**、测试和预测数据帧。**
事实上，测试数据集合是不是严格要求的**，因为H2O允许多重交叉验证自动执行。**
然而，正如上一段所暗示的，为了评估和比较装配模型的样品内性能**，我正在从列车数据中“切割”一个测试集。**
train_tbl <-
model_data_tbl %>%
filter(order_date <= "2007-03-19") test_tbl <-
model_data_tbl %>%
filter(order_date >= "2006-10-02" &
order_date <= "2007-03-19")train_tbl %>% head()
## # A tibble: 6 x 7
## order_date revenue trend trend_sqr rev_lag_13 rev_lag_52 season
## <date> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
## 1 2005-01-03 0 53 2809 0 0 0
## 2 2005-01-10 54013487\. 54 2916 45011429\. 58814754\. 1
## 3 2005-01-17 40984715\. 55 3025 30075259\. 13926869\. 1
## 4 2005-01-24 0 56 3136 0 0 0
## 5 2005-01-31 0 57 3249 0 0 0
## 6 2005-02-07 51927116\. 58 3364 51049952\. 55440318\. 1
创建预测数据集时的主要考虑围绕着对预测变量的可能值和水平进行计算假设****
- 说到
trend特性，我只是从 _model_data tbl 数据集中选择它们。它们是基于数字索引的，这样就很好了 - 假设有订单的周数几乎都是年复一年一致的(还记得探索性分析吗？)我将
season和rev_lag_52设置为一年前(52 周前)的值 rev_lag_13的值被设置为等于其上一季度的值(即 Q1 2007)
forecast_tbl <-
model_data_tbl %>%
filter(order_date > "2007-03-19") %>%
select(order_date:trend_sqr) %>%
cbind(season = model_data_tbl %>%
filter(between(order_date,
as.Date("2006-03-27"),
as.Date("2006-06-19"))) %>%
select(season),
rev_lag_52 = model_data_tbl %>%
filter(between(order_date,
as.Date("2006-03-27"),
as.Date("2006-06-19"))) %>%
select(rev_lag_52),
rev_lag_13 = model_data_tbl %>%
filter(between(order_date,
as.Date("2006-12-25"),
as.Date("2007-03-19"))) %>%
select(rev_lag_13)
) forecast_tbl %>% head()
## order_date revenue trend trend_sqr season rev_lag_52 rev_lag_13
## 1 2007-03-26 449709 169 28561 0 0.0 0
## 2 2007-04-02 0 170 28900 0 0.0 0
## 3 2007-04-09 89020602 171 29241 1 45948859.7 63603122
## 4 2007-04-16 70869888 172 29584 1 41664162.8 63305793
## 5 2007-04-23 0 173 29929 0 480138.8 0
## 6 2007-04-30 0 174 30276 0 0.0 0
和 H2O 一起做模特
终于准备好开始做模特了！
H2O是一个面向机器学习应用的高性能开源库，致力于分布式处理，这使得它适合于较小的内存项目，并且可以通过外部处理能力快速扩展以实现更大的目标。
它基于 Java，具有与 R 和 Python 的专用接口，并整合了许多监督和非监督的机器学习模型。在这个项目中，我特别关注 4 种算法:
- 广义线性模型(GLM)
- 梯度推进机(GBM)
- 我还使用了 AutoML 工具，并使用 leader 模型来比较性能
首先:启动一个H2O实例！
当 R 通过h2o.init命令启动 H2O 时，我可以指定内存分配池集群的大小。为了加快速度，我把它设置为“16G”。
h2o.init(max_mem_size = "16G")
## H2O is not running yet, starting it now...
##
## Note: In case of errors look at the following log files:
## C:\Users\LENOVO\AppData\Local\Temp\RtmpSWW88g/h2o_LENOVO_started_from_r.out
## C:\Users\LENOVO\AppData\Local\Temp\RtmpSWW88g/h2o_LENOVO_started_from_r.err
##
##
## Starting H2O JVM and connecting: . Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 4 seconds 712 milliseconds
## H2O cluster timezone: Europe/Berlin
## H2O data parsing timezone: UTC
## H2O cluster version: 3.26.0.10
## H2O cluster version age: 2 months and 4 days
## H2O cluster name: H2O_started_from_R_LENOVO_xwx278
## H2O cluster total nodes: 1
## H2O cluster total memory: 14.22 GB
## H2O cluster total cores: 4
## H2O cluster allowed cores: 4
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## H2O API Extensions: Amazon S3, Algos, AutoML, Core V3, TargetEncoder, Core V4
## R Version: R version 3.6.1 (2019-07-05)
我也喜欢关闭进度条，因为在某些情况下，输出消息可能会非常冗长。
h2o.no_progress()
下一步是安排响应和预测变量集。为了进行回归，您需要确保响应变量不是一个因子(否则H2O将进行分类)。
# response variable
y <- "revenue"
# predictors set: remove response variable and order_date from the set
x <- setdiff(names(train_tbl %>% as.h2o()), c(y, "order_date"))
随机森林
我将从安装一个random forest开始。
注意，我包括了nfolds参数。无论何时指定，该参数使交叉验证能够在不需要validation_frame的情况下执行——例如，如果设置为 5，它将执行 5 重交叉验证。
我还使用了一些控制参数来处理模型的运行时间:
- 我将
stopping_metric设置为RMSE作为提前停止的误差度量(当度量停止改善时，模型将停止构建新的树) - 通过
stopping_rounds，我指定了考虑提前停止前的训练轮数 - 我使用
stopping_tolerance来设置训练过程继续所需的最小改进
# random forest model
rft_model <-
h2o.randomForest(
x = x,
y = y,
training_frame = train_tbl %>% as.h2o(),
nfolds = 10,
ntrees = 500,
stopping_metric = "RMSE",
stopping_rounds = 10,
stopping_tolerance = 0.005,
seed = 1975
)
我现在用h2o.varimp_plot可视化变量重要性，它返回一个每个变量的排名贡献图，标准化为 0 到 1 之间的范围。
rft_model %>% h2o.varimp_plot()
model_summary功能允许访问关于模型参数的信息。
rft_model@model$model_summary
## Model Summary:
## number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1 27 27 12560 7
## max_depth mean_depth min_leaves max_leaves mean_leaves
## 1 14 10.29630 12 45 32.37037
这里我们可以看到，随机森林只使用了允许估计的最多 500 棵树中的 26 棵树(我用ntrees参数设置了这一点)。我们还可以估算出树的深度范围从 7 到 14(不是特别深的森林),每棵树的叶子数量范围从 12 到 45。
最后但同样重要的是，我可以用h2o.performance查看模型的性能
h2o.performance(rft_model, newdata = test_tbl %>% as.h2o())
## H2ORegressionMetrics: drf
##
## MSE: 3.434687e+13
## RMSE: 5860620
## MAE: 3635468
## RMSLE: 9.903415
## Mean Residual Deviance : 3.434687e+13
## R^2 : 0.9737282
该模型实现了 97.4% 的高R^2，这意味着特征变量的变化解释了响应变量的几乎所有可变性。
另一方面，RMSE显得相当大！RMSE 的高值可能是由于少量高误差预测的存在(如异常值的情况)，考虑到响应变量的波动性，这并不奇怪。
谈到基于误差的指标，如 RMSE、平均误差、均方误差等。没有好坏的绝对值，因为它们是用响应变量的单位表示的。通常，你想要实现一个更小的RMSE，因为这意味着更高的预测能力，但是对于这个项目，我将简单地使用这个指标来比较不同模型的相对性能。**
扩展到许多型号
让我们以编程的方式概括性能评估，一次性计算、评估和比较多个模型。
首先，我安装了几个型号，并确保我为所有型号启用了cross-validation。注意，对于 GBM ，我指定了与用于随机森林相同的参数，但是有大量的参数可以用来控制模型估计的几个方面(我不会涉及这些，因为这超出了本项目的范围)
# gradient boosting machine model
gbm_model <-
h2o.gbm(
x = x,
y = y,
training_frame = as.h2o(train_tbl),
nfolds = 10,
ntrees = 500,
stopping_metric = "RMSE",
stopping_rounds = 10,
stopping_tolerance = 0.005,
seed = 1975
)
# generalised linear model (a.k.a. elastic net model)
glm_model <-
h2o.glm(
x = x,
y = y,
training_frame = as.h2o(train_tbl),
nfolds = 10,
family = "gaussian",
seed = 1975
)
我还运行了非常方便的automl功能，可以安装多个模型并优化网格搜索。就像我对其他模型所做的那样，我可以指定一系列参数来指导这个函数，比如几个stopping度量和max_runtime_secs来节省计算时间。
automl_model <-
h2o.automl(
x = x,
y = y,
training_frame = as.h2o(train_tbl),
nfolds = 5,
stopping_metric = "RMSE",
stopping_rounds = 10,
stopping_tolerance = 0.005,
max_runtime_secs = 60,
seed = 1975
)
检查引导板将显示适合的模型
automl_model@leaderboard
## model_id
## 1 GBM_2_AutoML_20200112_111324
## 2 GBM_4_AutoML_20200112_111324
## 3 StackedEnsemble_BestOfFamily_AutoML_20200112_111324
## 4 GBM_1_AutoML_20200112_111324
## 5 DeepLearning_grid_1_AutoML_20200112_111324_model_1
## 6 DeepLearning_grid_1_AutoML_20200112_111324_model_5
##
## mean_residual_deviance rmse mse mae rmsle
## 1 1.047135e+14 10232963 1.047135e+14 5895463 NaN
## 2 1.070608e+14 10347021 1.070608e+14 5965855 NaN
## 3 1.080933e+14 10396794 1.080933e+14 5827000 NaN
## 4 1.102083e+14 10498016 1.102083e+14 5113982 NaN
## 5 1.104058e+14 10507419 1.104058e+14 6201525 NaN
## 6 1.146914e+14 10709407 1.146914e+14 6580217 NaN
##
## [22 rows x 6 columns]
如你所见，顶部的模型是一个梯度推进机模型。还有几个深度学习模型和一个堆叠合奏**，H2O迎战超级学习者。**
性能评价
首先，我将所有模型保存在一个文件夹中，这样我就可以访问它们，并通过一系列函数以编程方式处理性能指标
# set path to get around model path being different from project path
path = "/02_models/final/"
# Save GLM model
h2o.saveModel(glm_model, path)
# Save RF model
h2o.saveModel(rft_model, path)
# Save GBM model
h2o.saveModel(gbm_model, path)
# Extracs and save the leader autoML model
aml_model <- automl_model@leader
h2o.saveModel(aml_model, path)
可变重要性图
让我们从可变重要性图开始。之前我在random forest模型中使用了绘图功能，但是现在我想一次将它们全部绘制出来，这样我就可以比较和对比结果了。
有许多库(如 IML 、 PDP 、 VIP 、 DALEX 等更受欢迎的库)有助于机器学习模型可解释性**、特征解释和一般性能评估。在这个项目中，我使用的是vip包。**
这些库的主要优势之一是它们与其他 R 包如gridExtra的兼容性。
p_glm <- vip(glm_model) + ggtitle("GLM")
p_rft <- vip(rft_model) + ggtitle("RF")
p_gbm <- vip(gbm_model) + ggtitle("GBM")
p_aml <- vip(aml_model) + ggtitle("AML")
grid.arrange(p_glm, p_rft, p_gbm, p_aml, nrow = 2)
Seasonality和以前的收入水平(lags)在几乎所有车型中都排在前 3 位(唯一的例外是 GBM )。相反，没有一个模型发现trend及其对应的squared是反应变量变化的有力解释者。
性能指标
perf_gbm_model <-
h2o.performance(gbm_model, newdata = as.h2o(test_tbl))
perf_gbm_model
## H2ORegressionMetrics: gbm
##
## MSE: 1.629507e+13
## RMSE: 4036716
## MAE: 2150460
## RMSLE: 9.469847
## Mean Residual Deviance : 1.629507e+13
## R^2 : 0.9875359
然而，为了评估和比较模型的性能，我将重点关注RMSE和R^2^。
使用h20.metric函数可以一次提取所有性能指标，由于某种原因，这似乎不适用于H2ORegressionMetrics对象。
perf_gbm_model %>%
h2o.metric()
## Error in paste0("No ", metric, " for ",
## class(object)) : argument "metric" is missing, with
## no default
此外，在这种情况下，错误消息不是特别有用，因为“metric”参数是可选的，默认情况下应该返回所有指标。这个问题似乎与执行回归有关，因为它实际上可以很好地处理H2OClassificationMetrics对象。
幸运的是，提供了一些有用的助手函数来单独提取单个指标，它们工作正常！
perf_gbm_model %>% h2o.r2()
## [1] 0.9875359
perf_gbm_model %>% h2o.rmse()
## [1] 4036716
因此，我将使用这些单独的助手来编写一个小函数，该函数对测试数据上的所有模型运行预测，并返回一个包含所有性能指标的简便 tibble。****
performance_metrics_fct <- function(path, data_tbl) {
model_h2o <- h2o.loadModel(path)
perf_h2o <- h2o.performance(model_h2o, newdata = as.h2o(data_tbl))
R2 <- perf_h2o %>% h2o.r2()
RMSE <- perf_h2o %>% h2o.rmse()
tibble(R2, RMSE)
}
现在我可以将这个公式传递给purrr包中的map函数，以迭代计算并编译所有模型中的RMSE和R^2^。为了正确地识别每个模型，我还确保从路径中提取模型的名称。
perf_metrics_test_tbl <- fs::dir_info(path = "/02_models/final_models/") %>%
select(path) %>%
mutate(metrics = map(path, performance_metrics_fct, data_tbl = test_tbl),
path = str_split(path, pattern = "/", simplify = T)[,2]
%>% substr(1,3)) %>%
rename(model = path) %>%
unnest(cols = c(metrics))perf_metrics_test_tbl %>%
arrange(desc(R2))
model R2 RMSE
AML 0.9933358 2951704
GBM 0.9881890 3929538
DRF 0.9751434 5700579
GLM -0.0391253 36858064
所有基于树的模型都达到非常高的R^2，其中有 autoML 模型(这是一个 GBM，记得吗？)达到惊人的 99.3% 并达到最低RMSE。另一方面， GLM 得到一个负 R^2 。
负 R 并非闻所未闻:R 将模型的拟合与水平直线的拟合进行比较，并计算模型解释的方差与直线(零假设)解释的方差的比例。如果拟合实际上比仅仅拟合一条水平线差，那么 R 平方可以是负的。
实际与预测图
最后但同样重要的是，为了提供模型性能的额外和更直观的显示，我将绘制所有模型的实际与预测**。**
我正在使用一个类似于我用来计算性能指标的函数，因为基本原理是相同的。
predict_fct <- function(path, data_tbl) {
model_h2o <- h2o.loadModel(path)
pred_h2o <- h2o.predict(model_h2o, newdata = as.h2o(data_tbl))
pred_h2o %>%
as_tibble() %>%
cbind(data_tbl %>% select(order_date))
}
正如我之前所做的，我将公式传递给一个map函数来迭代计算，并使用所有模型的test数据子集来编译prediction。
validation_tmp <- fs::dir_info(path = "/02_models/final_models/") %>%
select(path) %>%
mutate(pred = map(path, predict_fct, data_tbl = test_tbl),
path = str_split(path, pattern = "/", simplify = T)[,2] %>%
substr(1,3)) %>%
rename(model = path)
然而，得到的validation_tmp是一个嵌套的 tibble，每个单元格中的预测都存储为列表。
validation_tmp
## # A tibble: 4 x 2
## model pred
## <chr> <list>
## 1 AML <df[,2] [25 × 2]>
## 2 DRF <df[,2] [25 × 2]>
## 3 GBM <df[,2] [25 × 2]>
## 4 GLM <df[,2] [25 × 2]>
这需要几个额外的操作来获得一个可用于绘图的形状:取消列表嵌套，围绕order_date旋转预测，并将收入添加为actual。
validation_tbl <-
validation_tmp %>%
unnest(cols = c(pred)) %>%
pivot_wider(names_from = model,
values_from = predict) %>%
cbind(test_tbl %>%
select(actual = revenue)) %>%
rename(date = order_date)
现在，我要把这个绘图函数直接写在里 plotly
validation_tbl %>%
plot_ly() %>%
add_lines(x = ~ date, y = ~ actual, name = 'Actual') %>%
add_lines(x = ~ date, y = ~ DRF, name = 'Random Forest',
line = list(dash = 'dot')) %>%
add_lines(x = ~ date, y = ~ GBM, name = 'Gradient Boosting Machine',
line = list(dash = 'dash')) %>%
add_lines(x = ~ date, y = ~ AML, name = 'Auto ML',
line = list(dash = 'dot')) %>%
add_lines(x = ~ date, y = ~ GLM, name = 'Generalised Linear Model',
line = list(dash = 'dash')) %>%
layout(title = 'Total Weekly Sales - Actual versus Predicted (various models)',
yaxis = list(title = 'Millions of Dollars'),
xaxis = list(title = ''),
legend = list(orientation = 'h')
)
除了 GLM 模型，它产生了一条看似平坦的预测线(还记得负面的 R ？)，所有模型都很好地捕捉到了系列中的波峰和波谷。预测仅开始错过最后 2 个尖峰周围的响应变化的全部范围。
预测
无需编写任何新函数，因为performance_metrics_fct和predict_fct也可用于预测。
首先，我看一下性能指标
perf_metrics_cast_tbl <- fs::dir_info(path = "/02_models/final_models/") %>%
select(path) %>%
mutate(metrics = map(path, performance_metrics_fct,
data_tbl = forecast_tbl),
path = str_split(path, pattern = "/", simplify = T)[,2]
%>% substr(1,3)) %>%
rename(model = path) %>%
unnest(cols = c(metrics)) perf_metrics_cast_tbl %>%
arrange(desc(R2))
model R2 RMSE
GBM 0.8678649 14544327
AML 0.8363565 16185792
DRF 0.8042526 17702414
GLM -0.0617160 41227664
有趣的是，在顶部位置有一点互换，与 autoML 模型相比,“手动”GBM 在预测中表现更好。与验证指标相比，所有模型的性能指标都恶化了。
然后，我计算预测…
cast_tbl <- fs::dir_info(path = "/02_models/final_models/") %>%
select(path) %>%
mutate(pred = map(path, predict_fct, data_tbl = forecast_tbl),
path = str_split(path, pattern = "/", simplify = T)[,2] %>%
substr(1,3)) %>%
rename(model = path) %>%
unnest(cols = c(pred)) %>%
pivot_wider(names_from = model, values_from = predict) %>%
cbind(forecast_tbl %>% select(actual = revenue)) %>%
rename(date = order_date)
…并想象它
cast_tbl %>%
plot_ly() %>%
add_lines(x = ~ date, y = ~ actual, name = 'Actual') %>%
add_lines(x = ~ date, y = ~ DRF, name = 'Random Forest',
line = list(dash = 'dot')) %>%
add_lines(x = ~ date, y = ~ GBM, name = 'Gradient Boosting Machine',
line = list(dash = 'dash')) %>%
add_lines(x = ~ date, y = ~ AML, name = 'Auto ML',
line = list(dash = 'dot')) %>%
add_lines(x = ~ date, y = ~ GLM, name = 'Generalised Linear Model',
line = list(dash = 'dash')) %>%
layout(title = 'Total Weekly Sales - Actual versus Forecast (various models)',
yaxis = list(title = 'Millions of Dollars'),
xaxis = list(title = ''),
legend = list(orientation = 'h')
)
除了不出所料产生持平预测的 GLM 之外，所有车型都继续很好地捕捉 2 周开工、2 周停工的购买模式。此外，所有模型预测都未能捕捉到所有 3 个峰值的反应变量运动的全部范围，这表明在 2007 年可能有另一种动态在起作用，而当前的预测者无法控制这种动态。
forecast_tbl %>%
select(-trend, -trend_sqr) %>%
tail(10)
order_date revenue season rev_lag_52 rev_lag_13
4 2007-04-16 70869888 1 41664162.8 63305793
5 2007-04-23 0 0 480138.8 0
6 2007-04-30 0 0 0.0 0
7 2007-05-07 78585882 1 41617508.0 64787291
8 2007-05-14 78797822 1 48403283.4 59552955
9 2007-05-21 0 1 0.0 0
10 2007-05-28 0 0 0.0 0
11 2007-06-04 0 0 45696327.2 0
12 2007-06-11 75486199 1 53596289.9 70430702
13 2007-06-18 86509530 1 774190.1 60094495
最后一件事:完成后不要忘记关闭H2O实例！
h2o.shutdown(prompt = FALSE)
结束语
在这个项目中，我经历了建立时间序列机器学习管道和生成每周收入预测所需的各个步骤。
特别是，我用TSstudio和进行了一个更“传统”的探索性时间序列分析，并利用我收集的洞察力创建了许多预测器**。然后我用开源库H2O训练并验证了一系列机器学习模型，并且使用性能指标和实际与预测图比较了模型的准确性。**
结论
这只是制作每周收入预测的第一次尝试，显然还有很大的改进空间。尽管如此，一旦你有了这样的建模和预测管道，创建和测试几个模型和不同的预测集就会变得更加容易和快速。
数据序列是人工生成的这一事实并不理想，因为它不一定包含您在现实生活数据集中会遇到的动态。尽管如此，这促使我发挥创造力，让整个练习变得更加愉快。
预测总收入可能不是最好的策略，例如，通过product line或country分解响应变量可能会带来更好、更准确的预测。这超出了本项目的范围，但也可能是未来项目的主题！
我从这个练习中学到的一件事是H2O是绝对聪明的**！它的设置快速、直观，并具有广泛的定制选项。 AutoML 非常棒，支持 R 、 Python 和 Java 并且任何人都可以免费使用它的事实给了像 Google AutoML 和 AWS SageMaker AutoML 这样的平台一个机会！**
代码库
完整的 R 代码可以在 我的 GitHub 简介 中找到
参考
- 对于 H2O 网站T5 H2O 网站
- 对于 H2O 文档 H2O 文档
- 有关在 R 中使用 R 进行时间序列分析和预测的详细讨论
- 关于 TSstudio 的介绍 [关于 TSstudio 包的介绍](https://diegousai.io/2019/12/time-series-machine-learning-analysis-and-demand-forecasting/Introduction for the TSstudio Package)
- T21 简介机器学习可解释性
原载于 2019 年 12 月 11 日https://diegousei . io。**
扫一扫
2284
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
