Keras RNN с ячейками LSTM для прогнозирования нескольких выходных временных рядов на основе множественных временных рядов

Я хотел бы моделировать RNN с ячейками LSTM, чтобы предсказать несколько временных рядов выходных данных на основе нескольких временных рядов ввода. Чтобы быть конкретным, у меня есть 4 выходных временных ряда, y1 [t], y2 [t], y3 [t], y4 [t], каждая из которых имеет длину 3000 (t = 0,..., 2999). У меня также есть 3 временных интервала ввода, x1 [t], x2 [t], x3 [t], и каждая из них имеет длину 3000 секунд (t = 0,..., 2999). Цель состоит в том, чтобы предсказать y1 [t],.. y4 [t], используя все входные временные ряды до этой текущей точки времени i.e.:

Для модели с долговременной памятью я создал модель RNN с сохранением состояния, следуя следующему. keras-stateful-lstme. Основное различие между моим случаем и keras-stateful-lstme заключается в том, что у меня есть:

Мой код запущен. Однако результат предсказания модели плох даже при простых данных. Поэтому я хотел бы спросить вас, не ошибаюсь ли я.

В примере с игрушками наши входные временные ряды представляют собой простые косинусные и знаковые волны:

Обратите внимание, что y1 в момент времени t является просто значением x1 при t - 2. Также обратите внимание, что y3 в момент времени t является просто значением x1 на предыдущем шаге.

Используя эти функции, я создал 100 наборов временных рядов y1, y2, y3, x1, x2, x3, x4. Половина из них идет на тренировочные данные, а оставшаяся половина переходит к тестовым данным.

Пример временного ряда y1,.. y4 и x1,..., x3 показаны ниже:

Создать модель с сохранением состояния RNN с 10 скрытыми нейронами LSTM

Теперь для обучения и проверки модели RNN используется следующий код:

Как вы можете видеть, потеря обучения НЕ уменьшается!

Поскольку целевые временные ряды 1 и 3 имеют очень простые отношения с входными временными рядами (y1 [t] = x1 [t-2], y3 [t] = x4 [t-3]), я ожидал бы идеального производительность прогнозирования. Однако тестирование R2 в каждую эпоху показывает, что это не так. R2 в конечную эпоху составляет всего около 0,2 и 0,36. Очевидно, что алгоритм не сходится. Я очень озадачен этим результатом. Пожалуйста, дайте мне знать, что мне не хватает, и почему алгоритм не сходится.

Ответы

Ответ 1

Начальная нота. Если временные ряды были короткими (например, T = 30), нам не требовалось бы сдержанное LSTM, и классический LSTM работал бы хорошо. В вопросе OP длина временных рядов равна T = 3000, поэтому обучение может быть очень медленным с классическим LSTM. Обучение будет улучшено путем сокращения временных рядов на части и использования LSTM с сохранением состояния.

Режим состояния с N = batch_size. Модифицированные модели сложны с Keras, потому что вам нужно быть осторожными, как вы сокращаете временные ряды и выбираете размер партии. В вопросе ОП размер выборки равен N = 100. Поскольку мы можем принять для обучения модели с партией из ста (это не большое число), мы выберем batch_size = 100.

Выбор batch_size = N упрощает обучение, потому что вам не нужно указывать reset внутри эпох (поэтому нет необходимости писать обратный вызов on_batch_begin).

Остается вопрос о сокращении временных рядов. Резка немного техническая, поэтому я написал функцию обертки, которая работает во всех случаях.

def stateful_cut(arr, batch_size, T_after_cut):
    if len(arr.shape) != 3:
        # N: Independent sample size,
        # T: Time length,
        # m: Dimension
        print("ERROR: please format arr as a (N, T, m) array.")

    N = arr.shape[0]
    T = arr.shape[1]

    # We need T_after_cut * nb_cuts = T
    nb_cuts = int(T / T_after_cut)
    if nb_cuts * T_after_cut != T:
        print("ERROR: T_after_cut must divide T")

    # We need batch_size * nb_reset = N
    # If nb_reset = 1, we only reset after the whole epoch, so no need to reset
    nb_reset = int(N / batch_size)
    if nb_reset * batch_size != N:
        print("ERROR: batch_size must divide N")

    # Cutting (technical)
    cut1 = np.split(arr, nb_reset, axis=0)
    cut2 = [np.split(x, nb_cuts, axis=1) for x in cut1]
    cut3 = [np.concatenate(x) for x in cut2]
    cut4 = np.concatenate(cut3)
    return(cut4)

С этого момента легко подготовить модель. Поскольку пример OP очень прост, нам не нужна дополнительная предварительная обработка или регуляризация. Я описываю, как действовать шаг за шагом (поскольку нетерпеливый, полностью автономный код доступен в самом конце этого сообщения).

Сначала мы загружаем данные и изменяем их с помощью функции-оболочки.

import numpy as np
from keras.models import Sequential
from keras.layers import Dense, LSTM, TimeDistributed
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

##
# Data
##
N = X_train.shape[0] # size of samples
T = X_train.shape[1] # length of each time series
batch_size = N # number of time series considered together: batch_size | N
T_after_cut = 100 # length of each cut part of the time series: T_after_cut | T
dim_in = X_train.shape[2] # dimension of input time series
dim_out = y_train.shape[2] # dimension of output time series

inputs, outputs, inputs_test, outputs_test = \
  [stateful_cut(arr, batch_size, T_after_cut) for arr in \
  [X_train, y_train, X_test, y_test]]

Затем мы скомпилируем модель с 4 входами, 3 выходами и 1 скрытым слоем, содержащим 10 узлов.

##
# Model
##
nb_units = 10

model = Sequential()
model.add(LSTM(batch_input_shape=(batch_size, None, dim_in),
               return_sequences=True, units=nb_units, stateful=True))
model.add(TimeDistributed(Dense(activation='linear', units=dim_out)))
model.compile(loss = 'mse', optimizer = 'rmsprop')

Мы обучаем модель без состояний сброса. Мы можем сделать это только потому, что мы выбрали batch_size = N.

##
# Training
##
epochs = 100

nb_reset = int(N / batch_size)
if nb_reset > 1:
    print("ERROR: We need to reset states when batch_size < N")

# When nb_reset = 1, we do not need to reinitialize states
history = model.fit(inputs, outputs, epochs = epochs, 
                    batch_size = batch_size, shuffle=False,
                    validation_data=(inputs_test, outputs_test))

Мы получаем эволюцию тренировочных/тестовых потерь следующим образом:

Теперь мы определяем "мим-модель", которая является апатридом, но содержит наши веса с состоянием. [Почему так? Прогнозирование с помощью модели stateful через model.predict требует полной партии в Keras, но у нас может не быть полной партии для прогнозирования...]

## Mime model which is stateless but containing stateful weights
model_stateless = Sequential()
model_stateless.add(LSTM(input_shape=(None, dim_in),
               return_sequences=True, units=nb_units))
model_stateless.add(TimeDistributed(Dense(activation='linear', units=dim_out)))
model_stateless.compile(loss = 'mse', optimizer = 'rmsprop')
model_stateless.set_weights(model.get_weights())

Наконец, мы можем показать наши невероятные прогнозы в наших длинных временных рядах y1, y2 и y3 (синий для истинного выхода, оранжевый для прогнозируемых выходов):

Для y1:

Для y2:

Для y3:

Заключение: он работает почти идеально, если только для 2-3 первых дат, когда сериал непредсказуем по определению. Мы не наблюдаем никаких всплесков при переходе из одной партии для следующей партии.

Гораздо больше. Когда N велико, мы хотели бы выбрать batch_size | N с batch_size < N. Я написал полный код в https://github.com/ahstat/deep-learning/blob/master/rnn/4_lagging_and_stateful.py (часть C и D). Этот путь github также показывает эффективность классического LSTM для коротких временных рядов (часть A) и неэффективность для длинных временных рядов (часть B). Вероятно, я напишу сообщение в блоге об этом в следующем месяце.

Полный автономный код

################
# Code from OP #
################
import numpy as np
def random_sample(len_timeseries=3000):
    Nchoice = 600
    x1 = np.cos(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
    x2 = np.cos(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
    x3 = np.sin(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
    x4 = np.sin(np.arange(0,len_timeseries)/float(1.0 + np.random.choice(Nchoice)))
    y1 = np.random.random(len_timeseries)
    y2 = np.random.random(len_timeseries)
    y3 = np.random.random(len_timeseries)
    for t in range(3,len_timeseries):
        ## the output time series depend on input as follows: 
        y1[t] = x1[t-2] 
        y2[t] = x2[t-1]*x3[t-2]
        y3[t] = x4[t-3]
    y = np.array([y1,y2,y3]).T
    X = np.array([x1,x2,x3,x4]).T
    return y, X
def generate_data(Nsequence = 1000):
    X_train = []
    y_train = []
    for isequence in range(Nsequence):
        y, X = random_sample()
        X_train.append(X)
        y_train.append(y)
    return np.array(X_train),np.array(y_train)

Nsequence = 100
prop = 0.5
Ntrain = int(Nsequence*prop)
X, y = generate_data(Nsequence)
X_train = X[:Ntrain,:,:]
X_test  = X[Ntrain:,:,:]
y_train = y[:Ntrain,:,:]
y_test  = y[Ntrain:,:,:] 

#X.shape = (N sequence, length of time series, N input features)
#y.shape = (N sequence, length of time series, N targets)
print(X.shape, y.shape)
# (100, 3000, 4) (100, 3000, 3)

####################
# Cutting function #
####################
def stateful_cut(arr, batch_size, T_after_cut):
    if len(arr.shape) != 3:
        # N: Independent sample size,
        # T: Time length,
        # m: Dimension
        print("ERROR: please format arr as a (N, T, m) array.")

    N = arr.shape[0]
    T = arr.shape[1]

    # We need T_after_cut * nb_cuts = T
    nb_cuts = int(T / T_after_cut)
    if nb_cuts * T_after_cut != T:
        print("ERROR: T_after_cut must divide T")

    # We need batch_size * nb_reset = N
    # If nb_reset = 1, we only reset after the whole epoch, so no need to reset
    nb_reset = int(N / batch_size)
    if nb_reset * batch_size != N:
        print("ERROR: batch_size must divide N")

    # Cutting (technical)
    cut1 = np.split(arr, nb_reset, axis=0)
    cut2 = [np.split(x, nb_cuts, axis=1) for x in cut1]
    cut3 = [np.concatenate(x) for x in cut2]
    cut4 = np.concatenate(cut3)
    return(cut4)

#############
# Main code #
#############
from keras.models import Sequential
from keras.layers import Dense, LSTM, TimeDistributed
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

##
# Data
##
N = X_train.shape[0] # size of samples
T = X_train.shape[1] # length of each time series
batch_size = N # number of time series considered together: batch_size | N
T_after_cut = 100 # length of each cut part of the time series: T_after_cut | T
dim_in = X_train.shape[2] # dimension of input time series
dim_out = y_train.shape[2] # dimension of output time series

inputs, outputs, inputs_test, outputs_test = \
  [stateful_cut(arr, batch_size, T_after_cut) for arr in \
  [X_train, y_train, X_test, y_test]]

##
# Model
##
nb_units = 10

model = Sequential()
model.add(LSTM(batch_input_shape=(batch_size, None, dim_in),
               return_sequences=True, units=nb_units, stateful=True))
model.add(TimeDistributed(Dense(activation='linear', units=dim_out)))
model.compile(loss = 'mse', optimizer = 'rmsprop')

##
# Training
##
epochs = 100

nb_reset = int(N / batch_size)
if nb_reset > 1:
    print("ERROR: We need to reset states when batch_size < N")

# When nb_reset = 1, we do not need to reinitialize states
history = model.fit(inputs, outputs, epochs = epochs, 
                    batch_size = batch_size, shuffle=False,
                    validation_data=(inputs_test, outputs_test))

def plotting(history):
    plt.plot(history.history['loss'], color = "red")
    plt.plot(history.history['val_loss'], color = "blue")
    red_patch = mpatches.Patch(color='red', label='Training')
    blue_patch = mpatches.Patch(color='blue', label='Test')
    plt.legend(handles=[red_patch, blue_patch])
    plt.xlabel('Epochs')
    plt.ylabel('MSE loss')
    plt.show()

plt.figure(figsize=(10,8))
plotting(history) # Evolution of training/test loss

##
# Visual checking for a time series
##
## Mime model which is stateless but containing stateful weights
model_stateless = Sequential()
model_stateless.add(LSTM(input_shape=(None, dim_in),
               return_sequences=True, units=nb_units))
model_stateless.add(TimeDistributed(Dense(activation='linear', units=dim_out)))
model_stateless.compile(loss = 'mse', optimizer = 'rmsprop')
model_stateless.set_weights(model.get_weights())

## Prediction of a new set
i = 0 # time series selected (between 0 and N-1)
x = X_train[i]
y = y_train[i]
y_hat = model_stateless.predict(np.array([x]))[0]

for dim in range(3): # dim = 0 for y1 ; dim = 1 for y2 ; dim = 2 for y3.
    plt.figure(figsize=(10,8))
    plt.plot(range(T), y[:,dim])
    plt.plot(range(T), y_hat[:,dim])
    plt.show()

## Conclusion: works almost perfectly.