用Python的长短期记忆神经网络进行时间序列预测

长短期记忆递归神经网络具有学习长的观察序列的潜力。

这对于时间序列预测似乎非常不错，并且事实的确可能是这样的。

在本教程中，你将了解，如何对于一个一步单变量时序预测问题开发一个LSTM预测模型。

完成本教程后，您将知道：

如何为预测问题制定性能基准。
如何为一步时间序列预测设计一个强大的测试框架。
如何准备数据，开发和评估用于时间序列预测的LSTM递归神经网络。

让我们开始吧。

更新于2017年5月：修复了invert_scale（）函数中的错误，谢谢Max。

教程概述

这是一个很大的话题，我们的教程将会覆盖很多内容，快准备好吧！

本教程分为9个部分; 他们是：

洗发水销售额数据集
测试设置
持续性模型预测
LSTM数据准备
LSTM模型开发
LSTM预测
完整的LSTM例子
开发稳健的结果
教程扩展

Python环境

本教程假设您已经安装了Python SciPy环境。本教程可以使用Python 2或3。

您必须安装了Keras(2.0或更高版本）和TensorFlow或Theano其中一个的后端

本教程还假设您已经安装了scikit-learn，Pandas，与NumPy和Matplotlib库

如果你的环境需要帮助，请看这个帖子：

如何使用Anaconda设置Python环境进行机器学习和深度学习

洗发水销售数据集

该数据集描述了3年期间洗发剂的月销售额。

这些单位是一个销售计数，有36个观测值。原始数据集归功于Makridakis，Wheelwright和Hyndman（1998）。

您可以在这里下载和了解更多关于数据集的信息。

将数据集下载到当前工作目录，名称为“ shampoo-sales.csv ”。请注意，您可能需要删除DataMarket添加的页脚信息。

下面的示例加载数据并绘制被加载的数据集的图形。

# load and plot dataset
from pandas import read_csv
from pandas import datetime
from matplotlib import pyplot
# 加载数据集
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# 总结第一行
print(series.head())
# 绘制折线图
series.plot()
pyplot.show()

运行该示例将数据集加载为Pandas系列并输出前5行。

Month
1901-01-01 266.0
1901-02-01 145.9
1901-03-01 183.1
1901-04-01 119.3
1901-05-01 180.3
Name: Sales, dtype: float64

然后创建系列的折线图，显示出明显的增长趋势。

实验测试设置

我们将洗发水销售数据集分为两部分：训练集和测试集。

前两年的数据将用于训练数据集，剩余的一年数据将用于测试集。

例如：

# 将数据分为训练和测试两部分
X = series.values
train, test = X[0:-12], X[-12:]

我们将通过使用训练集来开发模型并在测试集上做出预测。

滚动预测，也称为前向模型验证，在这里将会被用到。

测试集的每个时间步都会被同一个预测模型预测一次。然后测试集中每一个实际值都会被取出给预测模型使用，即对下一个时间步做出预测。

例如：

# 前向验证
history = [x for x in train]
predictions = list()
for i in range(len(test)):
    # 作出预测...

这模拟了一个真实世界的情景，每个月都有新的洗发水销售额产生，并用于下个月的预测。

最后，将收集关于测试数据集的所有预测，并计算误差分数以总结模型的性能。计算将使用均方根误差（RMSE），因为它会惩罚较大的错误，并得出与预测数据相同单位的分数，即月度洗发水销售额。

例如：

from sklearn.metrics import mean_squared_erro
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)

持续性模型预测

持续性预测模型对于带有线性增长趋势的时间序列是一个好的基准预测。

持续性预测是使用前一时间步（t-1）的观测值预测当前时间步（t）的观测值。

我们可以通过从训练数据和历史积累的历史数据中获取最后一个观测数据，并用它预测当前的时间步长来实现这一点。

例如：

# make prediction
yhat = history[-1]

我们将把所有的预测结果累加到一个数组中，以便将它们直接与测试数据集进行比较。

下面列出了如何用持续性模型预测洗发水销售数据集的完整示例。

from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_erro
from math import sqrt
from matplotlib import pyplot
# 加载数据集
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# 将数据分为训练和测试两部分
X = series.values
train, test = X[0:-12], X[-12:]
# 前向验证
history = [x for x in train]
predictions = list()
for i in range(len(test)):
    # 作出预测
    predictions.append(history[-1])
    # 测试集中的观测值
    history.append(test[i])
# 报告性能
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)
# 绘制观测值和预测值对比的折线图
pyplot.plot(test)
pyplot.plot(predictions)
pyplot.show()

运行该示例将基于测试数据集输出预测模型的RMSE，等于约136.

RMSE: 136.761

还会创建测试数据集（蓝色）与预测值（橙色）的折线图，显示了文中持续性预测模型的效果。

有关时间序列预测的持续性模型的更多信息，请参阅此文章：

如何使用Python进行时间序列预测的基线预测

现在我们已经有了数据集的性能基准，我们可以开始为数据开发一个LSTM模型

LSTM数据准备

在我们能够将LSTM模型拟合到数据集之前，我们必须对数据进行转换。

本节分为三个步骤：

将时间序列转化为监督学习问题
将时间序列转换为平稳时序。
将观察结果转换成具体的比例。

将时间序列转化为监督学习

Keras中的LSTM模型假定您的数据分为输入（X）和输出（y）。

对于时间序列问题，我们可以通过使用上一个时间步（t-1）的观测值作为输入，以当前时间步（t）的观测值作为输出来实现。

我们可以使用Pandas中的shift（）函数来实现这个功能，这个功能会将一系列的所有值按指定的位数推下去。我们需要一个位置的移位，这将成为输入变量。时间序列就是输出变量。

然后，我们可以将这两个系列连接起来，创建一个DataFrame，以供监督学习。下压的系列将在顶部有一个新的没有值的位置。NaN（非数字）值将被用在这个位置上。我们将用0值代替这些NaN值，LSTM模型将不得不学习“该系列的起始”或“这里没有数据”这样的情况，因为并没有观察到销售量为0的月份。

下面的代码定义了一个名为timeseries_to_supervised（）的辅助函数。它需要一个原始时间序列数据的NumPy数组和一个移位序列的滞后或数来创建并用作输入。

# 把一个序列作为一个监督学习问题
def timeseries_to_supervised(data, lag=1):
    df = DataFrame(data)
    columns = [df.shift(i) for i in range(1, lag+1)]
    columns.append(df)
    df = concat(columns, axis=1)
    df.fillna(0, inplace=True)
    return df

我们可以用我们加载的洗发水销售数据集来测试这个功能，并将其转换为监督学习问题。


from pandas import read_csv
from pandas import datetime
from pandas import DataFrame
from pandas import concat

# 把一个序列作为一个监督学习问题
def timeseries_to_supervised(data, lag=1):
    df = DataFrame(data)
    columns = [df.shift(i) for i in range(1, lag+1)]
    columns.append(df)
    df = concat(columns, axis=1)
    df.fillna(0, inplace=True)
    return df

# load dataset
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# 转换为监督学习问题
X = series.values
supervised = timeseries_to_supervised(X, 1)
print(supervised.head())

运行该示例将输出新监督学习问题的前5行。

            0           0
0    0.000000  266.000000
1  266.000000  145.899994
2  145.899994  183.100006
3  183.100006  119.300003
4  119.300003  180.300003

有关将时间序列问题转化为监督学习问题的更多信息，请参阅以下文章：

时间序列预测作为监督学习

将时间序列转换为平稳的

洗发水销售数据集是不平稳的

这意味着数据中有一个依赖于时间的结构。具体而言，数据呈现出上升趋势。

稳定的数据更易于建模，很可能会导致更准定的预测。

趋势可以从观测值中删除，然后再加回到预测值，以便将预测返回到原始的比例尺，并计算可比较的误差分数。

消除趋势的标准方法是差分化数据。这是从前一个时间步（t-1）的观察值减去当前的观测值（t）。这消除了趋势，我们留下了一个差分化系列，或从一个时间步的观测值到下一个时间步观测值的变化。

我们可以使用pandas中的diff（）函数自动实现这一点。另外，我们可以得到更好的粒度控制，并写我们自己的功能来做到这一点，在这种情况下，它的灵活性是首选。

下面是一个称为difference（）的函数，用于差分化时间序列。请注意，系列中的第一个观察值会被忽略，因为之前没有能用于计算的观测值。

# 创建一个差分化序列
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return Series(diff)

我们也需要转换这个过程，以便把对差分化序列的预测恢复到它们的原始尺度。

下面的函数inverse_difference（）就能达到。

# 转换差分化序列的值
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]

我们可以通过区分整个系列来测试这些函数，然后将其返回到原始尺度，如下所示：

from pandas import read_csv
from pandas import datetime
from pandas import Series
 
# 创建一个差分化序列
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return Series(diff)
 
# 转换差分化序列的值
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]
 
# 加载数据集
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
print(series.head())
# 将数据转化为平稳序列
differenced = difference(series, 1)
print(differenced.head())
# 反向转换
inverted = list()
for i in range(len(differenced)):
    value = inverse_difference(series, differenced[i], len(series)-i)
    inverted.append(value)
inverted = Series(inverted)
print(inverted.head())

运行示例输出加载数据的前5行，然后输出差分化序列的前5行，最后输出反向转换序列的前5行。

请注意，原始数据集中的第一个观测值已从差分化数据中移除。除此之外，最后一组数据与预期的一样匹配。

Month
1901-01-01    266.0
1901-02-01    145.9
1901-03-01    183.1
1901-04-01    119.3
1901-05-01    180.3
 
Name: Sales, dtype: float64
0   -120.1
1     37.2
2    -63.8
3     61.0
4    -11.8
dtype: float64
 
0    145.9
1    183.1
2    119.3
3    180.3
4    168.5
dtype: float64

有关平稳化时间序列和差分化序列的更多信息，请参阅文章：

将时间序列按比例缩放

像其他神经网络一样，LSTM希望数据大小能控制在神经网络使用的激活函数的范围内。

LSTM的默认激活函数是双曲正切（tanh），它输出-1和1之间的值。这是时间序列数据的首选范围。

为了使实验公平化，必须在训练数据集上计算缩放系数（最小值和最大值），并将其应用于缩放测试数据集和任何预测。这是为了避免使用来自测试数据集的信息影响实验，这样可能给模型带来一个小优势。

我们可以使用MinMaxScaler类将数据集转换为范围[-1，1] 。像其他scikit-learn转换类一样，它需要以行和列的矩阵格式提供数据。因此，我们必须在转换之前重塑我们的NumPy数组。

例如：

# 转换范围
X = series.values
X = X.reshape(len(X), 1)
scaler = MinMaxScaler(feature_range=(-1, 1))
scaler = scaler.fit(X)
scaled_X = scaler.transform(X)

再次，我们必须将预测的比例倒置，使得数值回到原始的比例，以便可以解释结果并计算出可比较的误差分数。

# 反向转换
inverted_X = scaler.inverse_transform(scaled_X)

把所有这些合并，下面的例子转换了洗发水销售数据的规模。

from pandas import read_csv
from pandas import datetime
from pandas import Series
from sklearn.preprocessing import MinMaxScale
# 加载数据
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
print(series.head())
# 按比例转换
X = series.values
X = X.reshape(len(X), 1)
scaler = MinMaxScaler(feature_range=(-1, 1))
scaler = scaler.fit(X)
scaled_X = scaler.transform(X)
scaled_series = Series(scaled_X[:, 0])
print(scaled_series.head())
# 反响转换
inverted_X = scaler.inverse_transform(scaled_X)
inverted_series = Series(inverted_X[:, 0])
print(inverted_series.head())

运行该是首先输出加载数据的前5行，然后输出缩放数据的前5行，然后反向缩放前5行，匹配原始数据。


Month
1901-01-01    266.0
1901-02-01    145.9
1901-03-01    183.1
1901-04-01    119.3
1901-05-01    180.3
 
Name: Sales, dtype: float64
0   -0.478585
1   -0.905456
2   -0.773236
3   -1.000000
4   -0.783188
dtype: float64
 
0    266.0
1    145.9
2    183.1
3    119.3
4    180.3
dtype: float

现在我们知道如何为LSTM网络准备数据，我们可以开始开发我们的模型。

LSTM模型开发

长短期记忆网络（LSTM）是一种递归神经网络（RNN）。

这种类型的网络的好处是，它可以学习和记忆长序列，并不依赖于预先指定的窗口滞后观察作为输入。

在Keras中，这被称为有状态，并且在定义LSTM层时涉及到将“ 有状态 ”参数设置为“ 真 ”。

默认情况下，Keras中的LSTM层在一个批处理数据之间保持状态。一组数据是训练数据集的固定大小的行数，它定义了在更新网络的权重之前需要处理多少模式。在默认情况下，批次之间的LSTM层的状态被清除，因此我们必须使LSTM有状态。通过调用reset_states()函数，这使我们可以细粒度地控制LSTM层的状态。

LSTM层期望输入在具有维度的矩阵中:[样本、时间步骤、特征].

样本：这些是来自域的独立观测值，通常是数据行。
时间步：对于给定的观察，这些是给定变量的单独的时间步。
特点：这是观察到的独立措施。

我们在如何为网络设计洗发水销售数据方面有一定的灵活性。我们将保持它的简单并构建出问题，因为时间序列中每一个时间步都是一个单独的样本，具有一个特征。

鉴于训练数据集被定义为x输入和y输出,它必须被重塑成样本/时间步/功能格式。例如:

X, y = train[:, 0:-1], train[:, -1]
X = X.reshape(X.shape[0], 1, X.shape[1])

输入数据的形状必须使用“batch_input_shape”参数作为tuple在LSTM层中指定，它指定要读取的每批数据的预期数量、时间步的数量和特性的数量。

批量大小通常比样本总数小很多。它和时间点的数量一起，定义了网络学习数据的速度(权重的更新频率)。

定义LSTM层的最后一个输入参数是神经元的数量，也称为内存单元或块的数量。这是一个相当简单的问题，1到5之间的数字就足够了。

下面的代码会创建一个单独的LSTM隐藏层，通过“ batch_input_shape ”参数指定输入层的期望值。

 layer = LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True)

该网络需要输出层中的单个神经元以线性激活来预测下一时间步洗发剂销量。

一旦指定了网络，就必须使用后端数学库（如TensorFlow或Theano）将其编译为高效的符号表示形式。

在编译网络时，我们必须指定一个损失函数和优化算法。我们将使用“ mean_squared_error ”作为损失函数，因为它与我们感兴趣的RMSE紧密匹配，以及高效的ADAM优化算法。

使用Sequential Keras API来定义网络，下面的代码片段将创建和编译网络。

model = Sequential()
model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')

一旦编译完毕，它可以与训练数据相匹配。因为网络是有状态的，所以当内部状态重置时，我们必须控制。因此，我们必须在每个时间步都手动管理训练过程。

默认情况下，一个时间点的样本在暴露在神经网络之前是会被搅乱的。同样，这对于LSTM来说是不可取的，因为我们希望网络通过观察的顺序来构建状态。我们可以通过将“shuffle”设置为“False”来禁用这些示例。

此外，默认情况下，神经网络在每个时间点末尾都会报告大量关于模型的学习进度和性能的调试信息。我们可以通过将“ verbose ”参数设置为“ 0 ” 的级别来禁用此功能。

然后，我们可以在训练时期结束时重置内部状态，为下一次训练迭代做好准备。

下面是一个将神经网络手动拟合到训练数据集的循环。

for i in range(nb_epoch):
    model.fit(X, y, epochs=1, batch_size=batch_size, verbose=0, shuffle=False)
    model.reset_states()

将所有这些组合在一起，我们可以定义一个名为fit_lstm()的函数，它将训练并返回一个LSTM模型。作为参数，它采用监督学习格式的训练数据集、批处理大小、多个时间点和一些神经元。

def fit_lstm(train, batch_size, nb_epoch, neurons):
    X, y = train[:, 0:-1], train[:, -1]
    X = X.reshape(X.shape[0], 1, X.shape[1])
    model = Sequential()
    model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
    model.add(Dense(1))
    model.compile(loss='mean_squared_error', optimizer='adam')
    for i in range(nb_epoch):
        model.fit(X, y, epochs=1, batch_size=batch_size, verbose=0, shuffle=False)
        model.reset_states()
    return model

batch_size必须设置为1.这是因为它必须是训练和测试数据集大小的一个因素。

模型的predict（）函数也受到批量大小的限制; 那么它必须设置为1，因为我们有兴趣对测试数据进行一步预测。

我们不会在本教程中调整网络参数。相反，我们将使用下面的配置，这是在一个小的尝试和错误中发现的：

批量大小：1
时间点数：3000
神经元：4

作为本教程的扩展，您可能希望探索不同的模型参数，并查看是否可以提高性能。

更新：请考虑尝试1500个时间点和1个神经元来预测，性能可能会更好！

接下来，我们将看看如何使用一个合适的LSTM模型来做出一个一步的预测。

LSTM预测

一旦LSTM模型拟合了训练数据，它就可以用来进行预测。

再次，我们有一些灵活性。我们可以决定在所有训练数据中一次拟合模型，然后根据测试数据每次预测一个新的时间步（我们称之为固定方法），或者我们可以重新每次拟合模型或更新模型将测试数据的步骤作为来自测试数据的新观测数据提供（我们称之为动态方法）。

在本教程中，我们将使用固定的方法来简化它，但是我们希望动态方法能够带来更好的模型技巧。

为了做出预测，我们可以在模型上调用predict（）函数。这需要3D NumPy数组输入作为参数。在这种情况下，它将是一个值的数组，在上一个时间步的观察。

predict（）函数返回预测的阵列，每个输入行对应一个。因为我们提供了一个单一的输入，输出将是一个2D NumPy数组。

我们可以在下面列出的名为forecast（）的函数中捕获这个行为。给定一个拟合模型，在拟合模型时使用的批量大小（例如1）和测试数据中的一行，函数将从测试行中分离出输入数据，对其进行重构，并将预测作为单个浮点值。

def forecast(model, batch_size, row):
    X = row[0:-1]
    X = X.reshape(1, 1, len(X))
    yhat = model.predict(X, batch_size=batch_size)
    return yhat[0,0]

在训练过程中，内部状态在每个时间点后重置。在预测的同时，我们不希望在过程中重置内部状态。事实上，我们希望模型能够在测试数据集的每个时间步中预测状态。

这提出了一个问题，即在预测测试数据集之前，什么样的神经网络是一个良好的初始状态。

在本教程中，我们通过对训练数据集中的所有样本进行预测来设定状态。在理论上，应该设置内部状态来预测下一个时间步。

我们现在拥有所有需要的东西来将一个LSTM网络模型拟合到洗发剂销量数据上并评估此模型的性能。

在下一节中，我们讲把这些所有需要的东西合并在一起。

完整的LSTM例子

在本节中，我们将一个LSTM网络模型拟合到洗发剂销量数据上并评估此模型。

这将涉及到前面各节的所有内容。内容很多，所以让我们回顾一下：

从CSV文件加载数据集。
转换数据集使其能够拟合LSTM模型，其中包括：
1. 将数据转化为监督学习问题。
2. 将数据转换成平稳的
3. 转换数据，使其具有从-1到1的比例。
将有状态的LSTM网络模型拟合到训练数据中。
在测试数据上评估静态LSTM模型。
报告模型预测的性能。

在这个例子中需要注意的事情：

为了简洁起见，将缩放和反缩放行为移到函数scale()和invert_scale()中。
测试数据是通过对训练数据的缩放比例进行缩放的，，以确保测试数据的最小/最大值不影响模型。
数据转换顺序的调整是因为方便起见，首先使数据平稳，接着监督学习问题，然后缩放。
我们可以很容易地收集观察结果，在我们进行的过程中进行验证，并将它们进行区分。

完整的例子如下所示。


from pandas import DataFrame
from pandas import Series
from pandas import concat
from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_erro
from sklearn.preprocessing import MinMaxScale
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from math import sqrt
from matplotlib import pyplot
import numpy
 
# 用于加载数据集的日期——时间解析函数
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')
 
# 把序列建构成一个监督学习问题
def timeseries_to_supervised(data, lag=1):
    df = DataFrame(data)
    columns = [df.shift(i) for i in range(1, lag+1)]
    columns.append(df)
    df = concat(columns, axis=1)
    df.fillna(0, inplace=True)
    return df
 
# 创建一个差分化序列
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return Series(diff)
 
# 反向转换差分化的值
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]
 
# 将测试集的数据缩放到[-1,1]范围内
def scale(train, test):
    # fit scale
    scaler = MinMaxScaler(feature_range=(-1, 1))
    scaler = scaler.fit(train)
    # transform train
    train = train.reshape(train.shape[0], train.shape[1])
    train_scaled = scaler.transform(train)
    # transform test
    test = test.reshape(test.shape[0], test.shape[1])
    test_scaled = scaler.transform(test)
    return scaler, train_scaled, test_scaled
 
# 反缩放预测值
def invert_scale(scaler, X, value):
    new_row = [x for x in X] + [value]
    array = numpy.array(new_row)
    array = array.reshape(1, len(array))
    inverted = scaler.inverse_transform(array)
    return inverted[0, -1]
 
# 将一个LSTM模型拟合到训练数据集上
def fit_lstm(train, batch_size, nb_epoch, neurons):
    X, y = train[:, 0:-1], train[:, -1]
    X = X.reshape(X.shape[0], 1, X.shape[1])
    model = Sequential()
    model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
    model.add(Dense(1))
    model.compile(loss='mean_squared_error', optimizer='adam')
    for i in range(nb_epoch):
        model.fit(X, y, epochs=1, batch_size=batch_size, verbose=0, shuffle=False)
        model.reset_states()
    return model
 
# 做出一步预测
def forecast_lstm(model, batch_size, X):
    X = X.reshape(1, 1, len(X))
    yhat = model.predict(X, batch_size=batch_size)
    return yhat[0,0]
 
# 加载数据集
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
 
# 把数据变得平稳
raw_values = series.values
diff_values = difference(raw_values, 1)
 
# 转换数据变成监督学习问题
supervised = timeseries_to_supervised(diff_values, 1)
supervised_values = supervised.values
 
# 把数据分成训练数据集和测试数据集
train, test = supervised_values[0:-12], supervised_values[-12:]
 
# 缩放数据
scaler, train_scaled, test_scaled = scale(train, test)
 
# 拟合模型
lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
# 预测训练集
train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
lstm_model.predict(train_reshaped, batch_size=1)
 
# 在测试数据集上的前向验证
predictions = list()
for i in range(len(test_scaled)):
    # 做出一步预测
    X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
    yhat = forecast_lstm(lstm_model, 1, X)
    # 反向缩放
    yhat = invert_scale(scaler, X, yhat)
    # 反向转换差分化数据
    yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
    # 存储预测值
    predictions.append(yhat)
    expected = raw_values[len(train) + i + 1]
    print('Month=%d, Predicted=%f, Expected=%f' % (i+1, yhat, expected))
 
# 报告模型性能
rmse = sqrt(mean_squared_error(raw_values[-12:], predictions))
print('Test RMSE: %.3f' % rmse)
# 绘制观测值和预测值对比的折线图
pyplot.plot(raw_values[-12:])
pyplot.plot(predictions)
pyplot.show()

运行该示例将输出测试数据集中12个月的每个月的观测值和预测值。

该示例还输出所有预测的RMSE。该模型显示月洗发水销售额的均方根误差为71.721，优于RMSE达到136.761的持续性预测模型。

随机数用于播种LSTM，因此，您可能会在模型的单一运行中得到不同的结果。在下一节中我们将进一步讨论这个问题。

Month=1, Predicted=351.582196, Expected=339.700000
Month=2, Predicted=432.169667, Expected=440.400000
Month=3, Predicted=378.064505, Expected=315.900000
Month=4, Predicted=441.370077, Expected=439.300000
Month=5, Predicted=446.872627, Expected=401.300000
Month=6, Predicted=514.021244, Expected=437.400000
Month=7, Predicted=525.608903, Expected=575.500000
Month=8, Predicted=473.072365, Expected=407.600000
Month=9, Predicted=523.126979, Expected=682.000000
Month=10, Predicted=592.274106, Expected=475.300000
Month=11, Predicted=589.299863, Expected=581.300000
Month=12, Predicted=584.149152, Expected=646.900000
Test RMSE: 71.721

测试数据（蓝色）与预测值（橙色）的折线图也被创建，大致展现了模型的性能水平。

作为补充说明，您可以做一个快速的实验来建立您对测试工具的信任以及所有的转换和反向转换。

在前向验证中注释符合LSTM模型的线：

yhat = forecast_lstm(lstm_model, 1, X)

并将其替换为以下内容：

yhat = y

这应该产生一个具有完美性能的模型（例如，将观测结果作为输出预测的模型）。

结果应该如下，表明如果LSTM模型能够完美地预测该序列，则逆变换和误差计算将正确显示。

Month=1, Predicted=339.700000, Expected=339.700000
Month=2, Predicted=440.400000, Expected=440.400000
Month=3, Predicted=315.900000, Expected=315.900000
Month=4, Predicted=439.300000, Expected=439.300000
Month=5, Predicted=401.300000, Expected=401.300000
Month=6, Predicted=437.400000, Expected=437.400000
Month=7, Predicted=575.500000, Expected=575.500000
Month=8, Predicted=407.600000, Expected=407.600000
Month=9, Predicted=682.000000, Expected=682.000000
Month=10, Predicted=475.300000, Expected=475.300000
Month=11, Predicted=581.300000, Expected=581.300000
Month=12, Predicted=646.900000, Expected=646.900000
Test RMSE: 0.000

开发稳健结果

神经网络的一个困难是，它们给出了不同起始条件下的不同结果。

一种方法可能是修复Keras使用的随机数种子，以确保结果是可重现的。另一种方法是使用不同的实验设置来控制随机的初始条件。

有关机器学习中的更多随机性，请参阅文章：

感受机器学习中的随机性

我们可以多次重复上一节的实验，然后将平均RMSE作为一个指示，说明该配置将如何在平均水平上执行看不见的数据。

这通常称为多次重复或多次重启。

我们可以将模型拟合和前向验证包装在固定数量的重复循环中。每次迭代都可以记录运行的RMSE。然后我们可以总结RMSE分数的分布。

# 重复试验
repeats = 30
error_scores = list()
for r in range(repeats):
    # 拟合模型
    lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
    # 预测训练集并创建预测环境
    train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
    lstm_model.predict(train_reshaped, batch_size=1)
    # 在测试数据集上进行前向验证
    predictions = list()
    for i in range(len(test_scaled)):
        # 做出一步预测
        X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
        yhat = forecast_lstm(lstm_model, 1, X)
        # 反向转换比例
        yhat = invert_scale(scaler, X, yhat)
        # 反向差分化数据
        yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
        # 存储预测值
        predictions.append(yhat)
    # 报告性能
    rmse = sqrt(mean_squared_error(raw_values[-12:], predictions))
    print('%d) Test RMSE: %.3f' % (r+1, rmse))
    error_scores.append(rmse)

数据准备将和以前一样。

我们将使用30次重复，因为这足以提供RMSE分数的良好分布。

完整的例子如下所示。


from pandas import DataFrame
from pandas import Series
from pandas import concat
from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_erro
from sklearn.preprocessing import MinMaxScale
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from math import sqrt
from matplotlib import pyplot
import numpy
 
# 用于加载数据集的日期——时间解析函数
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')
 
# 把一个序列构建为监督学习问题
def timeseries_to_supervised(data, lag=1):
    df = DataFrame(data)
    columns = [df.shift(i) for i in range(1, lag+1)]
    columns.append(df)
    df = concat(columns, axis=1)
    df.fillna(0, inplace=True)
    return df
 
# 创建一个差分化序列
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return Series(diff)
 
# 反向转换差分化的值
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]
 
# 将训练集和测试集缩放到[-1,1]的范围中
def scale(train, test):
    # 拟合比例
    scaler = MinMaxScaler(feature_range=(-1, 1))
    scaler = scaler.fit(train)
    # 转换训练集
    train = train.reshape(train.shape[0], train.shape[1])
    train_scaled = scaler.transform(train)
    # 转换测试集
    test = test.reshape(test.shape[0], test.shape[1])
    test_scaled = scaler.transform(test)
    return scaler, train_scaled, test_scaled
 
# 对预测值进行反缩放
def invert_scale(scaler, X, value):
    new_row = [x for x in X] + [value]
    array = numpy.array(new_row)
    array = array.reshape(1, len(array))
    inverted = scaler.inverse_transform(array)
    return inverted[0, -1]
 
# 把一个LSTM模型拟合到训练数据集上
def fit_lstm(train, batch_size, nb_epoch, neurons):
    X, y = train[:, 0:-1], train[:, -1]
    X = X.reshape(X.shape[0], 1, X.shape[1])
    model = Sequential()
    model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
    model.add(Dense(1))
    model.compile(loss='mean_squared_error', optimizer='adam')
    for i in range(nb_epoch):
        model.fit(X, y, epochs=1, batch_size=batch_size, verbose=0, shuffle=False)
        model.reset_states()
    return model
 
# 做出一步预测
def forecast_lstm(model, batch_size, X):
    X = X.reshape(1, 1, len(X))
    yhat = model.predict(X, batch_size=batch_size)
    return yhat[0,0]
 
# 加载数据集
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
 
# 将数据转换成平稳型
raw_values = series.values
diff_values = difference(raw_values, 1)
 
# 将数据转换为监督学习问题
supervised = timeseries_to_supervised(diff_values, 1)
supervised_values = supervised.values
 
# 将数据分为训练集和测试集
train, test = supervised_values[0:-12], supervised_values[-12:]
 
# 转换数据比例
scaler, train_scaled, test_scaled = scale(train, test)
 
# 重复试验
repeats = 30
error_scores = list()
for r in range(repeats):
    # 拟合模型
    lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
    # 预测训练集，创建预测环境
    train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
    lstm_model.predict(train_reshaped, batch_size=1)
    # 在测试集上面进行前向验证
    predictions = list()
    for i in range(len(test_scaled)):
        # 做出一步预测
        X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
        yhat = forecast_lstm(lstm_model, 1, X)
        # 反向转换比例
        yhat = invert_scale(scaler, X, yhat)
        # 反向转换差分化
        yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
        # 存储预测值
        predictions.append(yhat)
    # 报告性能
    rmse = sqrt(mean_squared_error(raw_values[-12:], predictions))
    print('%d) Test RMSE: %.3f' % (r+1, rmse))
    error_scores.append(rmse)
 
# 总结结果
results = DataFrame()
results['rmse'] = error_scores
print(results.describe())
results.boxplot()
pyplot.show()

运行该示例将输出每个重复的RMSE分数。运行结束提供收集的RMSE分数的汇总统计。

我们可以看到，月洗发剂销量的平均和标准偏差RMSE分别是138.491905和46.313783。

这是一个非常有用的结果，因为它表明上面报告的结果可能是一个统计巧合。实验表明，该模型可能与持续性模型平均一样好（136.761）。

这表明，至少需要进一步的模型调整。


1) Test RMSE: 136.191
2) Test RMSE: 169.693
3) Test RMSE: 176.553
4) Test RMSE: 198.954
5) Test RMSE: 148.960
6) Test RMSE: 103.744
7) Test RMSE: 164.344
8) Test RMSE: 108.829
9) Test RMSE: 232.282
10) Test RMSE: 110.824
11) Test RMSE: 163.741
12) Test RMSE: 111.535
13) Test RMSE: 118.324
14) Test RMSE: 107.486
15) Test RMSE: 97.719
16) Test RMSE: 87.817
17) Test RMSE: 92.920
18) Test RMSE: 112.528
19) Test RMSE: 131.687
20) Test RMSE: 92.343
21) Test RMSE: 173.249
22) Test RMSE: 182.336
23) Test RMSE: 101.477
24) Test RMSE: 108.171
25) Test RMSE: 135.880
26) Test RMSE: 254.507
27) Test RMSE: 87.198
28) Test RMSE: 122.588
29) Test RMSE: 228.449
30) Test RMSE: 94.427
             rmse
count   30.000000
mean   138.491905
std     46.313783
min     87.198493
25%    104.679391
50%    120.456233
75%    168.356040
max    254.507272

从下面显示的分布创建一个箱线图。这捕获了数据的中间以及范围和异常值结果。

这是一个实验设置，可用于比较LSTM模型的一个配置或设置到另一个配置。

教程扩展

本教程有很多扩展，我们可以考虑。

也许你可以自己探索一些，并在下面的评论中发表你的发现。

多步骤预测。实验设置可以改变，以预测下一个n阶段的步骤，而不是下一个单一的时间步骤。这也将允许更大的批量大小和更快的模型训练。请注意，在本教程中，尽管有新的观察值，并作为输入变量使用，我们基本上执行了一种12个一步的连续预测，模型并没有更新。
调整LSTM模型。模型没有调整;相反，这个配置是通过一些快速的尝试和错误发现的。我相信，至少可以通过调整神经元的数量和训练周期的数量来获得更好的结果。我还认为在训练模型期间通过回调提前停止可能会有用。
种子状态实验。在预测所有的培训数据是否有益之前，目前尚不清楚是否应当将该系统用于预测。这在理论上似乎是个好主意，但这需要加以论证。另外，也许在预测之前播种模型的其他方法是有益的。
更新模型。该模型可以再前向验证的每个时间步中进行更新。需要进行实验来确定是否从头开始重新构建模型会更好，或者用更多训练集（包括新样本））上的数据来更新权重。
输入时间步骤。LSTM输入支持样本的多个时间步。需要进行实验来观察是否包括滞后观测。
输入滞后功能。滞后观察可以被包括作为输入特征。需要实验来观察包括滞后特征是否提供任何好处，与AR（k）线性模型不同。
输入错误系列。可以构造一个错误序列（来自持续性模型的预测误差）并用作附加的输入特征，与MA（k）线性模型不同。需要进行实验，看看这是否能带来任何好处。
学习非固定。LSTM网络可能能够了解数据的趋势并做出合理的预测。需要进行实验来观察LSTM是否能够学习和有效地预测数据中留下的时间相关结构，如趋势和季节性。
对比无状态。本教程中使用了有状态的LSTM。结果应该与无状态的LSTM配置进行比较。
统计学意义。多次重复实验协议可以被进一步扩展以包括统计学意义测试，来证明不同构型的RMSE结果之间的差异是否具有统计学意义。

概要

在本教程中，您了解了如何开发用于时间序列预测的LSTM模型。

具体来说，你了解到：

如何准备用于开发LSTM模型的时间序列数据。
如何开发时间序列预测的LSTM模型。
如何使用强大的测试工具评估LSTM模型。