[译]Time Series Forecasting with the Long Short-Term Memory Network in Python

最新推荐文章于 2023-10-13 02:48:46 发布

EastWR

最新推荐文章于 2023-10-13 02:48:46 发布

阅读量4k

点赞数 5

分类专栏： deep-learning 文章标签： deep-learning LSTM

deep-learning 专栏收录该内容

4 篇文章

订阅专栏

[译]Time Series Forecasting with the Long Short-Term Memory Network in Python

写在开头的话，这篇文章是 Dr. Jason Brownlee 于2017-08-07发表的。是一篇非常优质的 LSTM 阅读文章。在网上没有找到合适的翻译，于是就自己动手了，^_^ 如有错误，欢迎指正。
附上原文

————- 分隔符以下是正文：

长短时记忆循环神经网络（LSTM）在学习观测长序列具有很好的前景。
这种方法看上去与解决时间序列预测问题完美的匹配，而事实上，它可能真的如此。
在本教程里，你将了解如何开发一个 LSTM 预测模型，以解决单变量单步时间序列预测问题。

阅读完本教程后，你会掌握：

对于预测问题，如何制定一个性能评判的基准（baseline）；
如何为单步时间序列预测问题设计一套鲁棒的测试过程；
针对时间序列预测问题，如何预处理数据，开发并评估一个 LSTM 循环神经网络；

让我们开始吧

Update May/2017: Fixed bug in invert_scale() function, thanks Max.

教程概览：

这篇文章涵盖面广。请读者阅读期间系好安全带，学习的小车准备开飞（滑稽脸）。

本教程会分为以下9个部分，它们是：

洗发水销售数据；
创建测试环节；
持久模型预测（译：原文Persistence Model Forecast
它的作用是直接通过训练数据产生baseline）；
LSTM 数据准备；
LSTM 模型开发；
LSTM 进行预测；
完成 LSTM 样例；
使预测结果具有鲁棒性；
扩展内容；

Python 环境

本教程假设你已经安装了Python SciPy。你可以使用Python2 或 3 都是ok 的。

你必须已经安装了Keras（2.0 或更高版本），同时安装了TensorFlow或Theano其中之一作为后台。

本教程假设你已经安装了scikit-learn， Pandas，NumPy 和 Matplotlib。

如果在配置开发环境方面需要帮助，请参考下面这篇文章：

How to Setup a Pyhon Environment for Machine Learning and Deep Learning with Anaconda

洗发水销售数据

这个数据集收集了三年的洗发水月销售量数据。
数据的单位是销售数量，一共包含36份观测数据。原始的数据集引用自Makridakis, Wheelwright, and Hyndman(1998)

你可以在这里下载并了解更多有关这个数据集的内容

下载数据集到你的本地目录你会得到一个名为’shampoo-sales.csv’的文件。注意你可能需要删除被数据网站添加的一些页脚信息。

加载数据集的样例代码如下：

# load and plot dataset
from pandas import read_csv
from pandas import datetime
from matplotlib import pyplot
# 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)
# summarize first few rows
print(series.head())
# line plot
series.plot()
pyplot.show()

运行样例程序，按照Pandas的Series格式读取并打印数据：

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

一个描述序列的图片被创建，显示出明显的上升趋势。

这里写图片描述

创建测试环节

我们会把洗发水销售数据切分成两部分：训练数据和测试数据。

前两年的数据会被切分成训练数据，剩下的数据作为测试集使用。

例如：

# split data into train and test
X = series.values
train, test = X[0:-12], X[-12:]

使用训练集建立模型，使用测试集对模型进行预测。

我们采用一种 ‘滚动预测’ 的场景，也被成为 ‘前行校验模型’。
测试数据集一次仅前进一个时间步长。模型会对当前时间步进行预测，然后测试集中实际的观测值会在下一个时间步的预测中使用。

例如：

# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
    # make prediction...

这种设置模拟了现实世界中每个月得到新的洗发水销售数据，然后再将新的数据用于预测下个月的预测的场景。

最终，测试集上所有的数据会被预测，然后计算汇总error score用以衡量模型的好坏。采用RMSE（均方误差）来评测模型因为它会惩罚大误差。然后会产生每个月的销售量预测只。

例如：

from sklearn.metrics import mean_squared_error
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_error
from math import sqrt
from matplotlib import pyplot
# 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)
# split data into train and test
X = series.values
train, test = X[0:-12], X[-12:]
# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
    # make prediction
    predictions.append(history[-1])
    # observation
    history.append(test[i])
# report performance
rmse = sqrt(mean_squared_error(test, predictions))
print('RMSE: %.3f' % rmse)
# line plot of observed vs predicted
pyplot.plot(test)
pyplot.plot(predictions)
pyplot.show()

运行样例程序，并打印 36个月的洗发水销售数据测试结果的RMSE值得到：

RMSE: 136.761

测试集（蓝色）和预测值（橙色）的折线图被创建，显示了持续预测模型中的上下文关系。

这里写图片描述

更多有关实现序列的持续预测问题，可以看文章

How to Make Baseline Predictions for Time Series Forecasting with Python

现在我们已经得到数据集的一个基准性能，现在我们可以使用数据集开始开发 LSTM 模型了。

LSTM 数据准备

在我们使用数据集拟合 LSTM 模型之前，我们必须转换数据的格式。

这部分内容被分成3个步骤：

将时间序列转换成有监督学习的问题；
转换时间序列数据，使它们是静态的；
转换观测值，使它们缩放到指定的尺度范围内；

将时间序列转换为有监督学习过程

在Keras中的 LSTM 模型会假设你的数据被分割成 input(X)和output(y)的形式。

对于时间序列问题，我们可以使上一个时间步(t-1)的观测值作为输入input，当前时刻(t)的观测值作为输出；

想要实现这一点我们可以使用Pandas自带的shift()函数，它可以使基于时间序列的值在指定距离下向后退。（译者注：shift函数可以使时间序列类型的数据，在以时间为index时，其所对应的值向过去或者未来滑动。）我们要求value与时间的对应关系改变1个距离单位，这会作为输入变量。时间序列原始对应的值作为输出变量。

我们可以连接这两组时间序列数据，使它们作为用于有监督学习的完备DataFrame数据。被向后移动的时间序列数据会空出一个空间来（这个空间原本是t时刻的值，现在它退到t-1时刻了，所以t时刻本身就没有值了），此时这个位置上不会被赋予一个值。一个NaN的标志（而不是一个数字）会被填充到这个位置上。我们会用0来替换这个NaN，这也导致 LSTM 不得不区分当数据集中一个月份对应0数值时，表示的是“这个序列的起始位置”还是“此处观测数据为0”两种情况。

下面的代码定义一个名为 timeseries_to_supervised() 的辅助函数来实现这个功能。它需要一个NumPy数组作为原始数据输入，并且赋予参数 lag 指明后退的距离。然后以得到有监督学习中的‘input’数据。

# frame a sequence as a supervised learning problem
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

# frame a sequence as a supervised learning problem
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)
# transform to supervised learning
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

更多有关将时间序列问题转化为有监督学习过程的内容可以参考下面的文章：

Time Series Forecasting as Supervised Learning

换转时间序列数据为平稳的（stationary）

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

这意味着该数据具有一种结构，它依赖于时间。更明确一点，数据中具有一种上升的趋势。

平稳的数据更容易建模也更容易产生精准的预测。

我们可以将这种 “趋势” 从观测数据中剔除出去，在随后得到的预测结果中再把这种‘趋势’添加回去，以使的预测的结果恢复和原始数据集一样的规模，然后再计算一个类似的错误分数。

消除趋势的标准方法是差值化数据。那就是从t-1时间步得到的观测数据是当前时刻t的观测数据做差得到的。这消除了’趋势’，并且我们得到了不同的时间序列，或者得到了从一个时间步到下一个时间步的变化。

我们可以使用pandas中的 diff() 函数来自动化的实现这一点。可选择的是，我们可以得到更加细粒度的控制，编写我们自己的函数来实现这一点，对于本例来说这是更好的。

下面是一个称为 difference()的函数，它是用来计算差值序列的。注意，第一个值被跳过了，因为它没有先验观测值来计算差值。

# create a differenced 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)

我们也需要反转这个过程，为了使基于差值序列得到的预测结果能够返回到它们的原始数据形式。

下面这个函数，称为 inverse_difference(), 逆向操作这个过程，然后返回它们的原始值，如下：

from pandas import read_csv
from pandas import datetime
from pandas import Series

# create a differenced 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)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]

# 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)
print(series.head())
# transform to be stationary
differenced = difference(series, 1)
print(differenced.head())
# invert transform
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条差值逆操作的前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的范围内。这对于时间序列数据是一个理想的范围。

为了使预测结果客观，缩放系数（min和max）值必须基于训练集进行计算，并将它们应用在测试集中，以及任何一次的预测结果中。这是为了避免从测试集中获取了知识污染了实验结果，这可能会导致试验结果有一个小的提升。

我们可以使用 MinMaxScaler class 来缩放数据到区间[-1,1]。和scikt-learn的其他转换类一样，它要求提供的数据是行列矩阵形式的。因此，我们在转换之前，必须先 reshape 我们的 NumPy 数组。

例如：

# transform scale
X = series.values
X = X.reshape(len(X), 1)
scaler = MinMaxScaler(feature_range=(-1, 1))
scaler = scaler.fit(X)
scaled_X = scaler.transform(X)

再一次，我们必须逆向缩放预测结果恢复到原始比例大小。以使结果是可解释的，并且有一个可比较的错误评分进行计算。

# invert transform
inverted_X = scaler.inverse_transform(scaled_X)

把所有的这些整合到一起，下面的例子就是转换洗发水销售数据的比例：

from pandas import read_csv
from pandas import datetime
from pandas import Series
from sklearn.preprocessing import MinMaxScaler
# 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)
print(series.head())
# transform scale
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())
# invert transform
inverted_X = scaler.inverse_transform(scaled_X)
inverted_series = Series(inverted_X[:, 0])
print(inverted_series.head())

运行上面的程序，打印前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: float64

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

LSTM 模型开发

长短时记忆网络（LSTM）是一种循环神经网络（RNN）。

这种网络的优点是能够学习和记住长序列，并且这种能力并不依赖于一个长的延迟观测窗口作为输入。

在Keras，这种能力被称为状态（stateful），当我们定义一个 LSTM 层时设置参数 stateful 为 True 即可。

默认情况下，在 Keras 中 LSTM 层会保存一个 batch 内部数据之间的状态。一个 bacth 的数据是训练集中一组固定行数的数据，它定义了在更新网络权重之间需要处理多少样例。在 LSTM 层中，不同batch之间的状态默认情况下(计算完后)会被清除，因此我们必须让 LSTM 保有状态性。这使得当 LSTM 层的状态被清除时，我们可以通过调用 reset_states()函数来对网络进行细粒度的控制。

LSTM 层期待输入的数组维度：[samples, time step, features]

Samples：它和观测数据本身是无关的，一般是数据的行数。
Time step：对于给定的观测数据，它们是给定变量下的独立时间步。
Features：在观测时间下，它们是一些独立的指标；

关于如何给神经网络构建洗发水销售数据，我们有一些灵活的方法。我们会保持它的简单性，在构建数据时，原始数据中每个时间步下的数据是一个单独的样本，并且一个时间步对应一个特征。（这里原文两个time step不知道怎么翻译好了，We will keep it simple and frame the problem as each time step in the original sequence is one separate sample, with one timestep and one feature.）

给定训练集被定义为 X 输入和 Y 输出的形式，必须 reshape 数据格式变为： Samples/TimeSteps/Features，例如：

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

输入数据的格式必须在 LSTM 层被指定，使用’batch_input_shape’参数，我们会传递一个元组赋予这个参数，以指明每读取一个batch时，期待多少条samples，多少个tiem steps以及多少个features。

batch的大小一般远小于样本的数量。它和epoch一起，决定了网络以多快的速度来学习数据（参数被多久更新一次）。

定义LSTM层的最后一个重要参数是，神经元的数量，也称为记忆单元或记忆块的数量。这是一个合理简单的小问题，一般取之1-5之间就足够了。

下面创建一个单独的LSTM隐层，通过 ‘batch_input_shape’ 参数指定了期待的输入。

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

网络要求输出只有一个神经元，使用 linear 激活函数，在下一个时间步中来产生洗发水销售数据的预测值。

网络一旦被定义，就需要使用后台库用高效的符号表示，例如 TensorFlow 或 Theano。

在编译网络时，我们必须指定一个损失函数和优化算法。我们会采用’mean_squared_error’来作为损失函数因为它和我们所感兴趣的 RMSE时接近的。同时，我们采用高效的ADAM优化算法。

使用Keras的Sequential模型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')

一但编译完成，它便可以拟合训练数据。因为网络是有状态的，所以当内部状态发生改变时，我们必须加以控制。因此，在贯穿整个预定义训练轮数（epoch）中的一轮训练时，我们必须人为的管理训练过程。

默认情况下，在一轮训练内，样本集在输送至网络前会被提前打乱顺序（译：洗牌操作 shuffle）。再次指出，对于LSTM而言这并不是我们所希望的，因为我们希望随着网络在贯穿学习观测序列数据时能够建立一种状态信息。我们可以通过设置 ‘shuffle’ 参数为 False 来禁用对样本的洗牌操作。

同样是默认的设置，在每一轮训练结束后，网络会输出很多关于学习过程和模型性能的调试信息。我们可以通过设置 ‘verbose’ 为第 ‘0’ 档，来禁止输出。

此时我们可以在每个训练轮次（epoch）结束后，重置网络状态，以备下一轮次的训练。

下面是一个循环，可以人为的控制网络使其拟合训练数据。

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模型。至于参数选择，则会以有监督学习的格式从训练数据中获取，batch 大小， epoch（训练轮数）数量，神经元数量。

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() 函数同样受batch size的约束；将batch size设置为1是因为我们的目的在于，在测试集上作出单步预测。

在本教程中，我们不会对网络进行调参；取而代之的是，我们会做下面一些工作，来发现一些小的改进和错误。

Batch Size： 1
Epochs： 3000
Neurons： 4

作为对本教程的一个扩展，你可能会感兴趣尝试不同的参数然后观察你是否可以改进模型性能。

更新：考虑尝试1500轮(epoch)训练和1个神经元，性能可能会更好！

下面，我们将会了解如何通过拟合好的LSTM模型来给出单步的预测结果。

LSTM 预测

一旦模型拟合完训练数据后，它可以被用来做预测。

此外，我们有一些灵活的选择。我们可以决定让模型一次性拟合全部的训练数据，然后在测试集中每次预测一个时间步（我们称这种方法为固定方式）。或者我们也可以随着测试集中数据可以作为观测数据使用时，在每一个时间步中重新拟合或者更新模型（我们称这种方法为动态方式）。

在本教程中，我们采用固定方式因为它比较简单。尽管，我们更期待采用动态方式，以得到更好的模型。

要进行预测时，我们可以调用模型的 predict() 函数。它要求一个三维的 NumPy 数组作为输入。在本用例中，输入的数据是从前一个时间步中，观测到的包含一个值的数组。

predict()函数会返回一个预测结果的数组，每一个值对应输入数据的一行。因为我们提供的单值输入，输出的结果将会是有一个值的 2 维 NumPy 数组。（原文：Because we are providing a single input, the output will be a 2D NumPy array with one value.）

我们可以通过下面列出的 forecast() 函数来展示这一过程。给定一个拟合好的模型，采用训练时设定好的 batch size（e.g.1），再加上测试集中的一行数据，这个函数会将输入数据进行分割，reshape it，然后返回一个浮点型的单值作为预测结果。

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() 函数中，这样会更简洁；
测试集上的缩放比例采用和拟合训练数据集一样的缩放尺度，以确保测试集上的最小/最大值不会对模型产生影响；
出于便利性的考虑，数据转换的顺序为，首先是数据稳定化，然后是转换有监督学习问题，最后是缩放数据；
数据差值会在切分数据为 train 和 test 集之间完成。在应用向前校验（walk-forward validation）时，我们本可以很容易收集数据并对它们进行差值计算。但是为了可读性，我决定不这样做了。

完整的例子如下所示：

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_error
from sklearn.preprocessing import MinMaxScaler
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

# date-time parsing function for loading the dataset
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')

# frame a sequence as a supervised learning problem
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

# create a differenced 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)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]

# scale train and test data to [-1, 1]
def scale(train, test):
    # fit scaler
    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

# inverse scaling for a forecasted value
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]

# fit an LSTM network to training data
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

# make a one-step forecast
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]

# load dataset
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)

# transform data to be stationary
raw_values = series.values
diff_values = difference(raw_values, 1)

# transform data to be supervised learning
supervised = timeseries_to_supervised(diff_values, 1)
supervised_values = supervised.values

# split data into train and test-sets
train, test = supervised_values[0:-12], supervised_values[-12:]

# transform the scale of the data
scaler, train_scaled, test_scaled = scale(train, test)

# fit the model
lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
# forecast the entire training dataset to build up state for forecasting
train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
lstm_model.predict(train_reshaped, batch_size=1)

# walk-forward validation on the test data
predictions = list()
for i in range(len(test_scaled)):
    # make one-step forecast
    X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
    yhat = forecast_lstm(lstm_model, 1, X)
    # invert scaling
    yhat = invert_scale(scaler, X, yhat)
    # invert differencing
    yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
    # store forecast
    predictions.append(yhat)
    expected = raw_values[len(train) + i + 1]
    print('Month=%d, Predicted=%f, Expected=%f' % (i+1, yhat, expected))

# report performance
rmse = sqrt(mean_squared_error(raw_values[-12:], predictions))
print('Test RMSE: %.3f' % rmse)
# line plot of observed vs predicted
pyplot.plot(raw_values[-12:])
pyplot.plot(predictions)
pyplot.show()

运行上面的结果，会打印每个月的期待值和测试集中的预测结果。

样例也会输出所有预测结果的 RMSE（均方根误差）。模型每个月份的洗发水销售数据 RMSE 值为 71.721，这要好于固定模型获得的 136.761 的RMSE值。

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

这样会产生一个完美熟练的模型。（e.g. 预测结果和观测结果完全一样）

这个结果应该如下所示，如果 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确保结果可复现的一个方法是固定随机数种子。另一种方法是使用不同的实验设置，以控制随机的初始条件。

更多关于机器学习去随机化的过程，可以参考下面这篇文章：

Embrace Randomness in Machine Learning

我们可以重复上面的实验多次，以获得一个平均 RMSE值，以衡量平均的情况下，模型在未知数据上性能的期望。

这常被称为多次重复，或多重启动。

我们可以将模型拟合和向前校验封装在一个固定重复的循环中。每一轮循环得到的 RMSE 值可以被记录下来。我们可以总结 RMSE 值的分布。

# repeat experiment
repeats = 30
error_scores = list()
for r in range(repeats):
    # fit the model
    lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
    # forecast the entire training dataset to build up state for forecasting
    train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
    lstm_model.predict(train_reshaped, batch_size=1)
    # walk-forward validation on the test data
    predictions = list()
    for i in range(len(test_scaled)):
        # make one-step forecast
        X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
        yhat = forecast_lstm(lstm_model, 1, X)
        # invert scaling
        yhat = invert_scale(scaler, X, yhat)
        # invert differencing
        yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
        # store forecast
        predictions.append(yhat)
    # report performance
    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_error
from sklearn.preprocessing import MinMaxScaler
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

# date-time parsing function for loading the dataset
def parser(x):
    return datetime.strptime('190'+x, '%Y-%m')

# frame a sequence as a supervised learning problem
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

# create a differenced 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)

# invert differenced value
def inverse_difference(history, yhat, interval=1):
    return yhat + history[-interval]

# scale train and test data to [-1, 1]
def scale(train, test):
    # fit scaler
    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

# inverse scaling for a forecasted value
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]

# fit an LSTM network to training data
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

# make a one-step forecast
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]

# load dataset
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)

# transform data to be stationary
raw_values = series.values
diff_values = difference(raw_values, 1)

# transform data to be supervised learning
supervised = timeseries_to_supervised(diff_values, 1)
supervised_values = supervised.values

# split data into train and test-sets
train, test = supervised_values[0:-12], supervised_values[-12:]

# transform the scale of the data
scaler, train_scaled, test_scaled = scale(train, test)

# repeat experiment
repeats = 30
error_scores = list()
for r in range(repeats):
    # fit the model
    lstm_model = fit_lstm(train_scaled, 1, 3000, 4)
    # forecast the entire training dataset to build up state for forecasting
    train_reshaped = train_scaled[:, 0].reshape(len(train_scaled), 1, 1)
    lstm_model.predict(train_reshaped, batch_size=1)
    # walk-forward validation on the test data
    predictions = list()
    for i in range(len(test_scaled)):
        # make one-step forecast
        X, y = test_scaled[i, 0:-1], test_scaled[i, -1]
        yhat = forecast_lstm(lstm_model, 1, X)
        # invert scaling
        yhat = invert_scale(scaler, X, yhat)
        # invert differencing
        yhat = inverse_difference(raw_values, yhat, len(test_scaled)+1-i)
        # store forecast
        predictions.append(yhat)
    # report performance
    rmse = sqrt(mean_squared_error(raw_values[-12:], predictions))
    print('%d) Test RMSE: %.3f' % (r+1, rmse))
    error_scores.append(rmse)

# summarize results
results = DataFrame()
results['rmse'] = error_scores
print(results.describe())
results.boxplot()
pyplot.show()

运行上面的程序，会打印每一轮循环得到的 RMSE值。程序末尾会给出一个 RMSE 值的统计分析。

我们可以看到，洗发水销售数据的 RMSE 的平均值和标准差分别是 138.491905 和 46.313783。

这是一个非常有用的结果，因为它显示出之前得到的结果很可能只是一个统计波动值。

这个结果表明，模型的平均性能和持久化模型的结果差不多一样好，如果没有变得更糟的话。

这至少显示出，模型的进一步调优是必须的。

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

一个箱形图如下，它描述了数据的中间值，取值范围，和离群结果。
（译：最上面的黑线是最大值，最小面的黑线是最小值。矩形的长表示标准差到均值点的波动范围，绿色的线是排名50%的数据点）

这里写图片描述

教程扩展

关门这篇教程，有很多地方值得我们扩展考虑。

或许，你可以自己去探索这些内容，然后在评论中发表你的发现。

多步预测：实验可以设置为预测接下来一次预测 n-time 时间步的结果，而不是一次只给出 1 个时间步的结果。这可能就允许设置更大的 batch size 和更快的训练过程。注意，在本教程中我们假定不更新模型，仅采用了一种类型的12个单步预测。尽管我们有新的观测数据可以利用并能够加入到新的输入数据中来。
调优LSTM模型：该模型是没有调优处理过的。参数的设置大多是一种带有错误的快速配置。我相信通过至少调节神经元数量和训练轮数将能够得到更好的结果。同时我也认为，通过回调函数提前停止训练过程会是有帮助的。
随机种子状态实现 Seed State Experiments：通过对训练数据进行预测，以获得一个预先的初始化模型，是否是有利的尚不明确。尽管在理论上这似乎是一个好方法，但我们需要将其展示出来。当然，或许有其他更好的方法来初始模型预测先验会更佳的有效。
更新模型：在每次向前校验过程中，可以更新模型。我们需要证明，是否需要重新拟合模型，或者在包含了一些新的样本后用更大的训练轮数来更新权重。
输入时间步：LSTM 模型支持一个样本内输出多个时间步。实验需要证明，是否提供一些延迟的观测数据作为额外的时间步是有益的。
输入延迟特征：一些延迟的观测或许可以作为特征加入到输入数据中来。需要采用实验来看这一步骤是否是有必须要的。不像一个 AR(k)线性模型。
输入错误序列：我们可以构建一个错误序列（从持久模型中得到的预测错误）作为额外的特征输入。这和MA(k)线性模型是不同的，需要我们采用实验去验证这一点
学习非稳定数据 Learn Non-Stationary：LSTM 模型或许可以学习数据的变化趋势并给出更佳合理的预测。实验需要去证明，是否具有时间依赖的结构能够被 LSTM 模型学习并有效利用。类似变化趋势，季节性。
对比无状态 Contrast Stateless：本教程使用的有状态 LSTM 模型。我们可以用无状态的做一个对比分析。
统计重要性 Statistical Significance：可以进行多次重复试验，来看看不同参数配置下 RMSE 值是否会有差异，且是否具有统计结果。