基于长短期记忆神经网络LSTM的多步长时间序列预测

345 阅读22分钟

长短时记忆网络(LSTM)是一种能够学习和预测长序列的递归神经网络。LSTMs除了学习长序列外,还可以学习一次多步预测,这对于时间序列的预测非常有用。LSTMs的一个困难在于,它们可能难以配置,而且需要大量的准备工作才能获得适合学习的格式的数据。

在本教程中,您将了解如何使用Keras在Python中开发用于多步骤时间序列预测的LSTM。

完成本教程后,您将知道:

  • 如何为多步时间序列预测准备数据。
  • 如何建立多步时间序列预测的LSTM模型。
  • 如何评价一个多步骤的时间序列预测。

教程概述

本教程分为四个部分。它们是:

  • 洗发水的销售数据集
  • 数据准备和模型评估
  • 持久性模型
  • 多步LSTM

环境

本教程假设您已经安装了Python SciPy环境。您可以在这个示例中使用python2或python3。

本教程假设您已经安装了Keras v2.0或更高版本,并安装了TensorFlow或Theano后端

本教程还假设您已经安装了scikit-learn、panda、NumPy和Matplotlib。

Shampoo Sales Dataset

这个数据集描述了3年期间洗发水的月销售量。单位是销售计数,有36个观察值。原始数据集归功于Makridakis, Wheelwright和Hyndman(1998)。您可以在这里下载并了解有关数据集的更多信息。下面的示例加载并创建已加载数据集的图。

接下来,我们将研究实验中使用的模型配置和测试工具。

数据准备和模型评估

本节描述本教程中使用的数据准备和模型评估

数据分割

我们将洗发水销售数据集分为两部分:培训和测试集。培训数据集采用前两年的数据,测试集采用后一年的数据。使用训练数据集开发模型,并对测试数据集进行预测。过去12个月的观察结果如下:

"3-01",339.7
"3-02",440.4
"3-03",315.9
"3-04",439.3
"3-05",401.3
"3-06",437.4
"3-07",575.5
"3-08",407.6
"3-09",682.0
"3-10",475.3
"3-11",581.3
"3-12",646.9

多步预测

我们将设计一个多步骤的预测。对于数据集最后12个月的某一个月,我们需要做3个月的预测。这是已知的历史观测(t-1, t-2,…t-n)预测t, t+1和t+2。具体来说,从第二年的12月开始,我们必须预测1月、2月和3月。从一月开始,我们必须预测二月、三月和四月。一直到10月,11月,12月的预测从第3年的9月。总共需要10个3个月的预测,如下:

Dec, Jan, Feb, Mar
Jan, Feb, Mar, Apr
Feb, Mar, Apr, May
Mar, Apr, May, Jun
Apr,  May, Jun, Jul
May, Jun, Jul, Aug
Jun, Jul, Aug, Sep
Jul, Aug, Sep, Oct
Aug, Sep, Oct, Nov
Sep, Oct, Nov, Dec

模型评价

将使用滚动预测场景,也称为前向模型验证。测试数据集的每个时间步骤都将一次执行一个。将使用一个模型对时间步骤进行预测,然后从测试集中获取下个月的实际期望值,并将其提供给模型,用于下一个时间步骤的预测。这模拟了一个真实的场景,在这个场景中,每个月都可以获得新的洗发水销售数据,并将其用于下一个月的预测。这将通过列车和测试数据集的结构进行模拟。将收集测试数据集上的所有预测,并计算错误得分,以总结模型对每个预测时间步骤的技能。使用均方根误差(RMSE)来惩罚较大的误差,得到的分数与预测数据的单位相同,即月度洗发水销售。

持久性模型

时间序列预测的一个好的基线是持久性模型。这是一个预测模型,在这个模型中,最后一个观测值被向前持久化。由于它的简单性,它通常被称为幼稚的预测。

准备数据

第一步是将数据从一个系列转换成一个有监督的学习问题。也就是从一组数字到一组输入和输出模式。我们可以使用一个预先准备好的函数series_to_supervised()来实现这一点。

# convert time series into supervised learning problem
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
 n_vars = 1 if type(data) is list else data.shape[1]
 df = DataFrame(data)
 cols, names = list(), list()
 # input sequence (t-n, ... t-1)
 for i in range(n_in, 0-1):
  cols.append(df.shift(i))
  names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
 # forecast sequence (t, t+1, ... t+n)
 for i in range(0, n_out):
  cols.append(df.shift(-i))
  if i == 0:
   names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
  else:
   names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
 # put it all together
 agg = concat(cols, axis=1)
 agg.columns = names
 # drop rows with NaN values
 if dropnan:
  agg.dropna(inplace=True)
 return agg

该函数可以通过传入加载的级数值,其中n_in值为1,n_out值为3来调用;例如:

supervised = series_to_supervised(raw_values, 13)

接下来,我们可以将监督学习数据集分为训练集和测试集。我们知道在这个表单中,最后10行包含最后一年的数据。这些行组成测试集,其余的数据组成训练数据集。我们可以将所有这些放在一个新函数中,该函数接受加载的系列和一些参数,并返回一个准备建模的训练集和测试集。

# transform series into train and test sets for supervised learning
def prepare_data(series, n_test, n_lag, n_seq):
 # extract raw values
 raw_values = series.values
 raw_values = raw_values.reshape(len(raw_values), 1)
 # transform into supervised learning problem X, y
 supervised = series_to_supervised(raw_values, n_lag, n_seq)
 supervised_values = supervised.values
 # split into train and test sets
 train, test = supervised_values[0:-n_test], supervised_values[-n_test:]
 return train, test

我们可以用洗发水数据集来测试。下面列出了完整的示例。

from pandas import DataFrame
from pandas import concat
from pandas import read_csv
from pandas import datetime
 
# date-time parsing function for loading the dataset
def parser(x):
 return datetime.strptime('190'+x, '%Y-%m')
 
# convert time series into supervised learning problem
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
 n_vars = 1 if type(data) is list else data.shape[1]
 df = DataFrame(data)
 cols, names = list(), list()
 # input sequence (t-n, ... t-1)
 for i in range(n_in, 0-1):
  cols.append(df.shift(i))
  names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
 # forecast sequence (t, t+1, ... t+n)
 for i in range(0, n_out):
  cols.append(df.shift(-i))
  if i == 0:
   names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
  else:
   names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
 # put it all together
 agg = concat(cols, axis=1)
 agg.columns = names
 # drop rows with NaN values
 if dropnan:
  agg.dropna(inplace=True)
 return agg
 
# transform series into train and test sets for supervised learning
def prepare_data(series, n_test, n_lag, n_seq):
 # extract raw values
 raw_values = series.values
 raw_values = raw_values.reshape(len(raw_values), 1)
 # transform into supervised learning problem X, y
 supervised = series_to_supervised(raw_values, n_lag, n_seq)
 supervised_values = supervised.values
 # split into train and test sets
 train, test = supervised_values[0:-n_test], supervised_values[-n_test:]
 return train, test
 
# load dataset
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# configure
n_lag = 1
n_seq = 3
n_test = 10
# prepare data
train, test = prepare_data(series, n_test, n_lag, n_seq)
print(test)
print('Train: %s, Test: %s' % (train.shape, test.shape))

运行该示例首先打印整个测试数据集,也就是最后10行。还打印了训练集、测试数据集的形状和大小。

[[ 342.3  339.7  440.4  315.9]
 [ 339.7  440.4  315.9  439.3]
 [ 440.4  315.9  439.3  401.3]
 [ 315.9  439.3  401.3  437.4]
 [ 439.3  401.3  437.4  575.5]
 [ 401.3  437.4  575.5  407.6]
 [ 437.4  575.5  407.6  682. ]
 [ 575.5  407.6  682.   475.3]
 [ 407.6  682.   475.3  581.3]
 [ 682.   475.3  581.3  646.9]]
Train: (234), Test: (104)

我们可以看到测试数据集第一行的单个输入值(第一列)与第二年12月份洗发水销售的观察结果相匹配:

"2-12",342.3

我们还可以看到,对于每个观察中的1个输入值和3个输出值,每行包含4列。

做预测

下一步是进行持久性预测。我们可以在一个名为persistence()的函数中轻松实现持久性预测,该函数执行最后一次观察和要持久化的预测步骤的数量。这个函数返回一个包含预测的数组。

make a persistence forecast
def persistence(last_ob, n_seq):
 return [last_ob for i in range(n_seq)]

然后我们可以为测试数据集中从第二年12月到第三年9月的每个时间步骤调用这个函数。下面是一个函数make_forecasts(),它执行此操作,并将数据集的训练、测试和配置作为参数,并返回预测列表。

# evaluate the persistence model
def make_forecasts(train, test, n_lag, n_seq):
 forecasts = list()
 for i in range(len(test)):
  X, y = test[i, 0:n_lag], test[i, n_lag:]
  # make forecast
  forecast = persistence(X[-1], n_seq)
  # store the forecast
  forecasts.append(forecast)
 return forecasts

我们可以这样调用这个函数:

forecasts = make_forecasts(train, test, 13)

评估预测

最后一步是评估这些预测。我们可以通过计算多步骤预测的每个时间步的RMSE来实现这一点,在本例中给出了3个RMSE得分。下面的函数evaluate_forecasts()计算并打印每个预测时间步骤的RMSE。

# evaluate the RMSE for each forecast time step
def evaluate_forecasts(test, forecasts, n_lag, n_seq):
 for i in range(n_seq):
  actual = test[:,(n_lag+i)]
  predicted = [forecast[i] for forecast in forecasts]
  rmse = sqrt(mean_squared_error(actual, predicted))
  print('t+%d RMSE: %f' % ((i+1), rmse))

我们可以这样调用这个函数:

evaluate_forecasts(test, forecasts, 13)

这也有助于在原始数据集的背景下绘制预测图,从而了解RMSE分数在背景下与问题的关系。我们可以先绘制整个洗发水数据集,然后将每个预测绘制为红线。下面的函数plot_forecasts()将创建并显示此图。

# plot the forecasts in the context of the original dataset
def plot_forecasts(series, forecasts, n_test):
 # plot the entire dataset in blue
 pyplot.plot(series.values)
 # plot the forecasts in red
 for i in range(len(forecasts)):
  off_s = len(series) - n_test + i
  off_e = off_s + len(forecasts[i])
  xaxis = [x for x in range(off_s, off_e)]
  pyplot.plot(xaxis, forecasts[i], color='red')
 # show the plot
 pyplot.show()

我们可以像下面这样调用这个函数。请注意,12个月的测试集中保留的观察值是12个,而上面使用的10个监督学习输入/输出模式的观察值是10个。

# plot forecasts
plot_forecasts(series, forecasts, 12)

通过将持久化预测与原始数据集中的实际持久化值连接起来,我们可以更好地绘制图表。这将需要将最后一个观测值添加到预测前面。下面是改进后的plot_forecasts()函数的更新版本。

# plot the forecasts in the context of the original dataset
def plot_forecasts(series, forecasts, n_test):
 # plot the entire dataset in blue
 pyplot.plot(series.values)
 # plot the forecasts in red
 for i in range(len(forecasts)):
  off_s = len(series) - 12 + i - 1
  off_e = off_s + len(forecasts[i]) + 1
  xaxis = [x for x in range(off_s, off_e)]
  yaxis = [series.values[off_s]] + forecasts[i]
  pyplot.plot(xaxis, yaxis, color='red')
 # show the plot
 pyplot.show()

完整的示例

我们可以把上述这些代码片段拼在一起,下面列出了多步骤持久性预测的完整代码示例:

from pandas import DataFrame
from pandas import concat
from pandas import read_csv
from pandas import datetime
from sklearn.metrics import mean_squared_error
from math import sqrt
from matplotlib import pyplot
 
# date-time parsing function for loading the dataset
def parser(x):
 return datetime.strptime('190'+x, '%Y-%m')
 
# convert time series into supervised learning problem
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
 n_vars = 1 if type(data) is list else data.shape[1]
 df = DataFrame(data)
 cols, names = list(), list()
 # input sequence (t-n, ... t-1)
 for i in range(n_in, 0-1):
  cols.append(df.shift(i))
  names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
 # forecast sequence (t, t+1, ... t+n)
 for i in range(0, n_out):
  cols.append(df.shift(-i))
  if i == 0:
   names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
  else:
   names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
 # put it all together
 agg = concat(cols, axis=1)
 agg.columns = names
 # drop rows with NaN values
 if dropnan:
  agg.dropna(inplace=True)
 return agg
 
# transform series into train and test sets for supervised learning
def prepare_data(series, n_test, n_lag, n_seq):
 # extract raw values
 raw_values = series.values
 raw_values = raw_values.reshape(len(raw_values), 1)
 # transform into supervised learning problem X, y
 supervised = series_to_supervised(raw_values, n_lag, n_seq)
 supervised_values = supervised.values
 # split into train and test sets
 train, test = supervised_values[0:-n_test], supervised_values[-n_test:]
 return train, test
 
# make a persistence forecast
def persistence(last_ob, n_seq):
 return [last_ob for i in range(n_seq)]
 
# evaluate the persistence model
def make_forecasts(train, test, n_lag, n_seq):
 forecasts = list()
 for i in range(len(test)):
  X, y = test[i, 0:n_lag], test[i, n_lag:]
  # make forecast
  forecast = persistence(X[-1], n_seq)
  # store the forecast
  forecasts.append(forecast)
 return forecasts
 
# evaluate the RMSE for each forecast time step
def evaluate_forecasts(test, forecasts, n_lag, n_seq):
 for i in range(n_seq):
  actual = test[:,(n_lag+i)]
  predicted = [forecast[i] for forecast in forecasts]
  rmse = sqrt(mean_squared_error(actual, predicted))
  print('t+%d RMSE: %f' % ((i+1), rmse))
 
# plot the forecasts in the context of the original dataset
def plot_forecasts(series, forecasts, n_test):
 # plot the entire dataset in blue
 pyplot.plot(series.values)
 # plot the forecasts in red
 for i in range(len(forecasts)):
  off_s = len(series) - n_test + i - 1
  off_e = off_s + len(forecasts[i]) + 1
  xaxis = [x for x in range(off_s, off_e)]
  yaxis = [series.values[off_s]] + forecasts[i]
  pyplot.plot(xaxis, yaxis, color='red')
 # show the plot
 pyplot.show()
 
# load dataset
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# configure
n_lag = 1
n_seq = 3
n_test = 10
# prepare data
train, test = prepare_data(series, n_test, n_lag, n_seq)
# make forecasts
forecasts = make_forecasts(train, test, n_lag, n_seq)
# evaluate forecasts
evaluate_forecasts(test, forecasts, n_lag, n_seq)
# plot forecasts
plot_forecasts(series, forecasts, n_test+2)

运行该示例首先为每个预测的时间步骤打印RMSE。这为我们提供了LSTM在每个时间步上的性能基线,我们希望LSTM能够表现得更好。

t+1 RMSE: 144.535304
t+2 RMSE: 86.479905
t+3 RMSE: 121.149168

同时还绘制了具有多步持续性预测的原始时间序列图。这些线连接到每个预测的适当输入值。这个上下文说明持久性预测实际上是多么幼稚。

多步LSTM网络

在本节中,我们将使用持久性示例作为起点,并研究将LSTM适合于培训数据并对测试数据集进行多步骤预测所需的更改。

准备数据

在我们使用这些数据来培训LSTM之前,必须准备好这些数据。具体来说,还需要两项改动:

平稳性:数据显示出增长的趋势,必须通过差分来消除。

归一化:数据的规模必须缩小到-1到1之间的值,即LSTM单元的激活函数。

我们可以引入一个函数使数据平稳,称为差分()。这将把一系列的值转换成一系列的差异,一个更简单的表示。

# 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)

我们可以使用sklearn库中的MinMaxScaler来缩放数据。将这些组合在一起,我们可以更新prepare_data()函数,使其首先对数据进行差异处理并对其进行重新缩放,然后将转换转换为有监督学习问题,并像以前在持久性示例中所做的那样培训测试集。该函数现在除了返回火车和测试数据集之外,还返回一个标量。

# transform series into train and test sets for supervised learning
def prepare_data(series, n_test, n_lag, n_seq):
 # extract raw values
 raw_values = series.values
 # transform data to be stationary
 diff_series = difference(raw_values, 1)
 diff_values = diff_series.values
 diff_values = diff_values.reshape(len(diff_values), 1)
 # rescale values to -11
 scaler = MinMaxScaler(feature_range=(-11))
 scaled_values = scaler.fit_transform(diff_values)
 scaled_values = scaled_values.reshape(len(scaled_values), 1)
 # transform into supervised learning problem X, y
 supervised = series_to_supervised(scaled_values, n_lag, n_seq)
 supervised_values = supervised.values
 # split into train and test sets
 train, test = supervised_values[0:-n_test], supervised_values[-n_test:]
 return scaler, train, test

我们可以这样调用:

# prepare data
scaler, train, test = prepare_data(series, n_test, n_lag, n_seq)

拟合LSTM网络

接下来,我们需要将LSTM网络模型与训练数据相匹配。这首先要求训练数据集从2D数组[样本,特征]转换为3D数组[样本,时间步长,特征]。我们将把时间步骤固定在1,所以这个更改很简单。接下来,我们需要设计一个LSTM网络。我们将使用一个简单的结构,一个隐藏层和一个LSTM单元,然后是一个线性激活的输出层和三个输出值。该网络将采用均方误差损失函数和高效的亚当优化算法。LSTM是有状态的;这意味着我们必须在每个训练阶段结束时手动重置网络状态。这个网络将适用于1500个时代。对于训练和预测,必须使用相同的批大小,并且我们要求在测试数据集的每个步骤中都进行预测。这意味着必须使用批大小为1的批处理。批量大小为1也称为在线学习,因为每次训练模式结束后,网络权重都会在训练过程中更新(而不是小批量或批量更新)。我们可以将所有这些放在一个名为fit_lstm()的函数中。该函数接受一些关键参数,这些参数可用于以后优化网络,并返回一个适合进行预测的LSTM模型。

# fit an LSTM network to training data
def fit_lstm(train, n_lag, n_seq, n_batch, nb_epoch, n_neurons):
 # reshape training into [samples, timesteps, features]
 X, y = train[:, 0:n_lag], train[:, n_lag:]
 X = X.reshape(X.shape[0], 1, X.shape[1])
 # design network
 model = Sequential()
 model.add(LSTM(n_neurons, batch_input_shape=(n_batch, X.shape[1], X.shape[2]), stateful=True))
 model.add(Dense(y.shape[1]))
 model.compile(loss='mean_squared_error', optimizer='adam')
 # fit network
 for i in range(nb_epoch):
  model.fit(X, y, epochs=1, batch_size=n_batch, verbose=0, shuffle=False)
  model.reset_states()
 return model

我们可以这样调用:

# fit model
model = fit_lstm(train, 13115001)

网络配置未调优;如果您喜欢,可以尝试不同的参数。

LSTM预测

下一步是利用fit LSTM网络进行预测。使用合适的LSTM网络,可以通过调用model.predict()进行单个预测。同样,必须将数据格式化为具有[示例、时间步骤、特性]格式的3D数组。我们可以将它封装到一个名为forecast_lstm()的函数中。

make one forecast with an LSTM,
def forecast_lstm(model, X, n_batch):
 # reshape input pattern to [samples, timesteps, features]
 X = X.reshape(11len(X))
 # make forecast
 forecast = model.predict(X, batch_size=n_batch)
 # convert to array
 return [x for x in forecast[0, :]]

我们可以从make_forecasts()函数调用这个函数,并将其更新为接受模型作为参数。更新后的版本如下:

# evaluate the persistence model
def make_forecasts(model, n_batch, train, test, n_lag, n_seq):
 forecasts = list()
 for i in range(len(test)):
  X, y = test[i, 0:n_lag], test[i, n_lag:]
  # make forecast
  forecast = forecast_lstm(model, X, n_batch)
  # store the forecast
  forecasts.append(forecast)
 return forecasts

更新版本可以这样调用:

make forecasts
forecasts = make_forecasts(model, 1, train, test, 13)

逆变换

在做出预测之后,我们需要对转换进行反向操作,以便将值返回到原始的比例。这是必要的,以便我们可以计算出与其他模型(如上面的持久性预测)相比较的错误得分和图表。我们可以使用提供inverse_transform()函数的MinMaxScaler对象来直接反转预测的规模。我们可以通过将最近的观察值(前几个月的洗发水销售)添加到第一个预测值,然后将该值向下传播,来逆转这种差异。这有点复杂;我们可以将行为包装在一个函数名inverse_difference()中,该函数接受预测和预测之前的最后一个观测值作为参数,并返回反向预测。

# invert differenced forecast
def inverse_difference(last_ob, forecast):
 # invert first forecast
 inverted = list()
 inverted.append(forecast[0] + last_ob)
 # propagate difference forecast using inverted first value
 for i in range(1len(forecast)):
  inverted.append(forecast[i] + inverted[i-1])
 return inverted

把这些放在一起,我们可以创建一个inverse_transform()函数,该函数在每次预测中工作,首先反转尺度,然后反转差异,将预测返回到它们原来的尺度。

# inverse data transform on forecasts
def inverse_transform(series, forecasts, scaler, n_test):
 inverted = list()
 for i in range(len(forecasts)):
  # create array from forecast
  forecast = array(forecasts[i])
  forecast = forecast.reshape(1len(forecast))
  # invert scaling
  inv_scale = scaler.inverse_transform(forecast)
  inv_scale = inv_scale[0, :]
  # invert differencing
  index = len(series) - n_test + i - 1
  last_ob = series.values[index]
  inv_diff = inverse_difference(last_ob, inv_scale)
  # store
  inverted.append(inv_diff)
 return inverted

可以这样调用:

# inverse transform forecasts and test
forecasts = inverse_transform(series, forecasts, scaler, n_test+2)

我们还可以对输出部分测试数据集进行逆变换,从而正确计算RMSE分数,如下图所示:

actual = [row[n_lag:] for row in test]
actual = inverse_transform(series, actual, scaler, n_test+2)

我们也可以简化RMSE分数的计算,期望测试数据只包含输出值,如下所示:

def evaluate_forecasts(test, forecasts, n_lag, n_seq):
 for i in range(n_seq):
  actual = [row[i] for row in test]
  predicted = [forecast[i] for forecast in forecasts]
  rmse = sqrt(mean_squared_error(actual, predicted))
  print('t+%d RMSE: %f' % ((i+1), rmse))

完整的示例

我们可以将所有这些部分连接在一起,并将LSTM网络用于多步时间序列预测问题,下面提供了完整的代码清单:

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
from numpy import array
 
# date-time parsing function for loading the dataset
def parser(x):
 return datetime.strptime('190'+x, '%Y-%m')
 
# convert time series into supervised learning problem
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
 n_vars = 1 if type(data) is list else data.shape[1]
 df = DataFrame(data)
 cols, names = list(), list()
 # input sequence (t-n, ... t-1)
 for i in range(n_in, 0-1):
  cols.append(df.shift(i))
  names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
 # forecast sequence (t, t+1, ... t+n)
 for i in range(0, n_out):
  cols.append(df.shift(-i))
  if i == 0:
   names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
  else:
   names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
 # put it all together
 agg = concat(cols, axis=1)
 agg.columns = names
 # drop rows with NaN values
 if dropnan:
  agg.dropna(inplace=True)
 return agg
 
# 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)
 
# transform series into train and test sets for supervised learning
def prepare_data(series, n_test, n_lag, n_seq):
 # extract raw values
 raw_values = series.values
 # transform data to be stationary
 diff_series = difference(raw_values, 1)
 diff_values = diff_series.values
 diff_values = diff_values.reshape(len(diff_values), 1)
 # rescale values to -11
 scaler = MinMaxScaler(feature_range=(-11))
 scaled_values = scaler.fit_transform(diff_values)
 scaled_values = scaled_values.reshape(len(scaled_values), 1)
 # transform into supervised learning problem X, y
 supervised = series_to_supervised(scaled_values, n_lag, n_seq)
 supervised_values = supervised.values
 # split into train and test sets
 train, test = supervised_values[0:-n_test], supervised_values[-n_test:]
 return scaler, train, test
 
# fit an LSTM network to training data
def fit_lstm(train, n_lag, n_seq, n_batch, nb_epoch, n_neurons):
 # reshape training into [samples, timesteps, features]
 X, y = train[:, 0:n_lag], train[:, n_lag:]
 X = X.reshape(X.shape[0], 1, X.shape[1])
 # design network
 model = Sequential()
 model.add(LSTM(n_neurons, batch_input_shape=(n_batch, X.shape[1], X.shape[2]), stateful=True))
 model.add(Dense(y.shape[1]))
 model.compile(loss='mean_squared_error', optimizer='adam')
 # fit network
 for i in range(nb_epoch):
  model.fit(X, y, epochs=1, batch_size=n_batch, verbose=0, shuffle=False)
  model.reset_states()
 return model
 
# make one forecast with an LSTM,
def forecast_lstm(model, X, n_batch):
 # reshape input pattern to [samples, timesteps, features]
 X = X.reshape(11len(X))
 # make forecast
 forecast = model.predict(X, batch_size=n_batch)
 # convert to array
 return [x for x in forecast[0, :]]
 
# evaluate the persistence model
def make_forecasts(model, n_batch, train, test, n_lag, n_seq):
 forecasts = list()
 for i in range(len(test)):
  X, y = test[i, 0:n_lag], test[i, n_lag:]
  # make forecast
  forecast = forecast_lstm(model, X, n_batch)
  # store the forecast
  forecasts.append(forecast)
 return forecasts
 
# invert differenced forecast
def inverse_difference(last_ob, forecast):
 # invert first forecast
 inverted = list()
 inverted.append(forecast[0] + last_ob)
 # propagate difference forecast using inverted first value
 for i in range(1len(forecast)):
  inverted.append(forecast[i] + inverted[i-1])
 return inverted
 
# inverse data transform on forecasts
def inverse_transform(series, forecasts, scaler, n_test):
 inverted = list()
 for i in range(len(forecasts)):
  # create array from forecast
  forecast = array(forecasts[i])
  forecast = forecast.reshape(1len(forecast))
  # invert scaling
  inv_scale = scaler.inverse_transform(forecast)
  inv_scale = inv_scale[0, :]
  # invert differencing
  index = len(series) - n_test + i - 1
  last_ob = series.values[index]
  inv_diff = inverse_difference(last_ob, inv_scale)
  # store
  inverted.append(inv_diff)
 return inverted
 
# evaluate the RMSE for each forecast time step
def evaluate_forecasts(test, forecasts, n_lag, n_seq):
 for i in range(n_seq):
  actual = [row[i] for row in test]
  predicted = [forecast[i] for forecast in forecasts]
  rmse = sqrt(mean_squared_error(actual, predicted))
  print('t+%d RMSE: %f' % ((i+1), rmse))
 
# plot the forecasts in the context of the original dataset
def plot_forecasts(series, forecasts, n_test):
 # plot the entire dataset in blue
 pyplot.plot(series.values)
 # plot the forecasts in red
 for i in range(len(forecasts)):
  off_s = len(series) - n_test + i - 1
  off_e = off_s + len(forecasts[i]) + 1
  xaxis = [x for x in range(off_s, off_e)]
  yaxis = [series.values[off_s]] + forecasts[i]
  pyplot.plot(xaxis, yaxis, color='red')
 # show the plot
 pyplot.show()
 
# load dataset
series = read_csv('shampoo-sales.csv', header=0, parse_dates=[0], index_col=0, squeeze=True, date_parser=parser)
# configure
n_lag = 1
n_seq = 3
n_test = 10
n_epochs = 1500
n_batch = 1
n_neurons = 1
# prepare data
scaler, train, test = prepare_data(series, n_test, n_lag, n_seq)
# fit model
model = fit_lstm(train, n_lag, n_seq, n_batch, n_epochs, n_neurons)
# make forecasts
forecasts = make_forecasts(model, n_batch, train, test, n_lag, n_seq)
# inverse transform forecasts and test
forecasts = inverse_transform(series, forecasts, scaler, n_test+2)
actual = [row[n_lag:] for row in test]
actual = inverse_transform(series, actual, scaler, n_test+2)
# evaluate forecasts
evaluate_forecasts(actual, forecasts, n_lag, n_seq)
# plot forecasts
plot_forecasts(series, forecasts, n_test+2)

运行该示例首先为每个预测的时间步骤打印RMSE。我们可以看到,每个预测时间步的得分都比持久性预测好,在某些情况下要好得多。这表明配置的LSTM确实具有解决问题的技能。有趣的是,RMSE并没有像预期的那样随着预测范围的延长而逐渐变差。这是因为t+2比t+1更容易预测。这可能是因为向下的勾号比本系列中提到的向上勾号更容易预测(可以通过对结果进行更深入的分析来证实这一点)。

t+1 RMSE: 95.973221
t+2 RMSE: 78.872348
t+3 RMSE: 105.613951

还创建了一个带有预测(红色)的系列线图(蓝色)。从图中可以看出,虽然模型的技巧比较好,但是一些预测并不十分好,还有很大的改进空间。

扩展

  • 更新LSTM。在提供新数据时,更改示例以重新安装或更新LSTM。10个训练阶段应该足以用新的观察来重新训练。
  • 优化LSTM。网格搜索本教程中使用的一些LSTM参数,例如纪元数、神经元数和层数,看看是否可以进一步提高性能。
  • Seq2Seq。使用LSTMs的编码器-解码器范式预测每个序列,看看这是否有任何好处。
  • 时间范围。尝试预测不同的时间范围,看看网络的行为在不同的交货时间是如何变化的。

作者:沂水寒城,CSDN博客专家,个人研究方向:机器学习、深度学习、NLP、CV

Blog: yishuihancheng.blog.csdn.net

赞 赏 作 者

Python中文社区作为一个去中心化的全球技术社区,以成为全球20万Python中文开发者的精神部落为愿景,目前覆盖各大主流媒体和协作平台,与阿里、腾讯、百度、微软、亚马逊、开源中国、CSDN等业界知名公司和技术社区建立了广泛的联系,拥有来自十多个国家和地区数万名登记会员,会员来自以工信部、清华大学、北京大学、北京邮电大学、中国人民银行、中科院、中金、华为、BAT、谷歌、微软等为代表的政府机关、科研单位、金融机构以及海内外知名公司,全平台近20万开发者关注。

推荐阅读:

一文读懂高并发情况下的常见缓存问题\

用 Django 开发基于以太坊智能合约的 DApp\

一文读懂 Python 分布式任务队列 celery\

5 分钟解读 Python 中的链式调用\

用 Python 创建一个比特币价格预警应用

▼点击成为社区会员   喜欢就点个在看吧