LSTM时间序列预测周围的预测区间

数据挖掘 机器学习 深度学习 时间序列 预言 lstm
2021-10-08 00:55:35

有没有一种方法可以计算来自 LSTM(或其他递归)神经网络的时间序列预测的预测区间(概率分布)?

例如,假设我预测未来 10 个样本(t+1 到 t+10),基于最后 10 个观察到的样本(t-9 到 t),我预计 t+1 的预测会更多比 t+10 时的预测准确。通常,人们可能会在预测周围绘制误差线以显示间隔。使用 ARIMA 模型(假设误差为正态分布),我可以计算每个预测值周围的预测区间(例如 95%)。我可以从 LSTM 模型计算相同的(或与预测间隔相关的)吗?

我一直在使用 Keras/Python 中的 LSTM,遵循来自 machinelearningmastery.com 的大量示例我的示例代码(如下)基于这些示例。我正在考虑将问题重新定义为离散箱的分类,因为这会产生每个类的置信度,但这似乎是一个糟糕的解决方案。

有几个类似的主题(如下所示),但似乎没有任何东西可以直接解决 LSTM(或实际上其他)神经网络的预测间隔问题:

https://stats.stackexchange.com/questions/25055/how-to-calculate-the-confidence-interval-for-time-series-prediction

使用 ARIMA 与 LSTM 进行时间序列预测

from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from math import sin
from matplotlib import pyplot
import numpy as np

# Build an LSTM network and train
def fit_lstm(X, y, batch_size, nb_epoch, neurons):
    X = X.reshape(X.shape[0], 1, X.shape[1]) # add in another dimension to the X data
    y = y.reshape(y.shape[0], y.shape[1])      # but don't add it to the y, as Dense has to be 1d?
    model = Sequential()
    model.add(LSTM(neurons, batch_input_shape=(batch_size, X.shape[1], X.shape[2]), stateful=True))
    model.add(Dense(y.shape[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=1, shuffle=False)
        model.reset_states()
    return model

# Configuration
n = 5000    # total size of dataset
SLIDING_WINDOW_LENGTH = 30
SLIDING_WINDOW_STEP_SIZE = 1
batch_size = 10
test_size = 0.1 # fraction of dataset to hold back for testing
nb_epochs = 100 # for training
neurons = 8 # LSTM layer complexity

# create dataset
#raw_values = [sin(i/2) for i in range(n)]  # simple sine wave
raw_values = [sin(i/2)+sin(i/6)+sin(i/36)+np.random.uniform(-1,1) for i in range(n)]  # double sine with noise
#raw_values = [(i%4) for i in range(n)] # saw tooth

all_data = np.array(raw_values).reshape(-1,1) # make into array, add anothe dimension for sci-kit compatibility

# data is segmented using a sliding window mechanism
all_data_windowed = [np.transpose(all_data[idx:idx+SLIDING_WINDOW_LENGTH]) for idx in np.arange(0,len(all_data)-SLIDING_WINDOW_LENGTH, SLIDING_WINDOW_STEP_SIZE)]
all_data_windowed = np.concatenate(all_data_windowed, axis=0).astype(np.float32)

# split data into train and test-sets
# round datasets down to a multiple of the batch size
test_length = int(round((len(all_data_windowed) * test_size) / batch_size) * batch_size)
train, test = all_data_windowed[:-test_length,:], all_data_windowed[-test_length:,:]
train_length = int(np.floor(train.shape[0] / batch_size)*batch_size) 
train = train[:train_length,...]

half_size = int(SLIDING_WINDOW_LENGTH/2) # split the examples half-half, to forecast the second half
X_train, y_train = train[:,:half_size], train[:,half_size:]
X_test, y_test = test[:,:half_size], test[:,half_size:]

# fit the model
lstm_model = fit_lstm(X_train, y_train, batch_size=batch_size, nb_epoch=nb_epochs, neurons=neurons)

# forecast the entire training dataset to build up state for forecasting
X_train_reshaped = X_train.reshape(X_train.shape[0], 1, X_train.shape[1])
lstm_model.predict(X_train_reshaped, batch_size=batch_size)

# predict from test dataset
X_test_reshaped = X_test.reshape(X_test.shape[0], 1, X_test.shape[1])
yhat = lstm_model.predict(X_test_reshaped, batch_size=batch_size)

#%% Plot prediction vs actual

x_axis_input = range(half_size)
x_axis_output = [x_axis_input[-1]] + list(half_size+np.array(range(half_size)))

fig = pyplot.figure()
ax = fig.add_subplot(111)
line1, = ax.plot(x_axis_input,np.zeros_like(x_axis_input), 'r-')
line2, = ax.plot(x_axis_output,np.zeros_like(x_axis_output), 'o-')
line3, = ax.plot(x_axis_output,np.zeros_like(x_axis_output), 'g-')
ax.set_xlim(np.min(x_axis_input),np.max(x_axis_output))
ax.set_ylim(-4,4)
pyplot.legend(('Input','Actual','Predicted'),loc='upper left')
pyplot.show()

# update plot in a loop
for idx in range(y_test.shape[0]):

    sample_input = X_test[idx]
    sample_truth = [sample_input[-1]] + list(y_test[idx]) # join lists
    sample_predicted = [sample_input[-1]] + list(yhat[idx])

    line1.set_ydata(sample_input)
    line2.set_ydata(sample_truth)
    line3.set_ydata(sample_predicted)
    fig.canvas.draw()
    fig.canvas.flush_events()

    pyplot.pause(.25)
4个回答

直接来说,这是不可能的。但是,如果您以不同的方式对其进行建模,则可以得出置信区间。您可以代替正态回归将其视为估计连续概率分布。通过对每一步执行此操作,您可以绘制分布图。这样做的方法是内核混合网络(https://janvdvegt.github.io/2017/06/07/Kernel-Mixture-Networks.html,公开,我的博客)或密度混合网络(http://www.cedar .buffalo.edu/~srihari/CSE574/Chap5/Chap5.7-MixDensityNetworks.pdf),第一个使用内核作为基础并估计这些内核的混合,第二个估计分布的混合,包括每个分布的参数分布。您使用对数似然来训练模型。

对不确定性进行建模的另一种选择是在训练期间和推理期间使用 dropout。您会多次执行此操作,并且每次从您的后部获取样本时。你不会得到分布,只有样本,但它是最容易实现的,而且效果很好。

在您的情况下,您必须考虑生成 t+2 到 t+10 的方式。根据您当前的设置,您可能必须从上一个时间步进行采样并将其提供给下一个时间步。这对于第一种方法和第二种方法都不是很好。如果每个时间步有 10 个输出(t+1 到 t+10),那么所有这些方法都更干净,但不太直观。

保形预测作为一个流行词可能对您来说很有趣,因为它可以在许多条件下工作 - 特别是它不需要正态分布误差并且它几乎适用于任何机器学习模型。

Scott LocklinHenrik Linusson给出了两个很好的介绍

我将有点分歧,并认为计算置信区间在实践中通常不是一件有价值的事情。原因是你总是需要做出一大堆假设。即使对于最简单的线性回归,您也需要

  • 线性关系。
  • 多元正态性。
  • 没有或很少有多重共线性。
  • 没有自相关。
  • 同方差性。

更实用的方法是进行蒙特卡罗模拟。如果您已经知道或愿意对输入变量的分布做出假设,请获取一大堆样本并将其提供给您的 LSTM,现在您可以凭经验计算您的“置信区间”。

是的你可以。您唯一需要更改的是损失函数。实现分位数回归中使用的损失函数并将其集成。此外,您还想看看如何评估这些间隔。为此,我会使用 ICP、MIL 和 RMIL 指标。