Interwał prognozy wokół prognozy szeregów czasowych LSTM

14

Czy istnieje metoda obliczania przedziału predykcji (rozkładu prawdopodobieństwa) wokół prognozy szeregów czasowych z sieci neuronowej LSTM (lub innej cyklicznej)?

Powiedzmy na przykład, że przewiduję 10 próbek w przyszłości (t + 1 do t + 10), w oparciu o 10 ostatnio zaobserwowanych próbek (t-9 do t), oczekiwałbym, że przewidywanie przy t + 1 będzie większe dokładne niż prognoza przy t + 10. Zazwyczaj można narysować słupki błędów wokół prognozy, aby pokazać interwał. Za pomocą modelu ARIMA (przy założeniu błędów normalnie rozłożonych) mogę obliczyć przedział predykcji (np. 95%) wokół każdej przewidywanej wartości. Czy mogę obliczyć to samo (lub coś, co dotyczy przedziału prognoz) na podstawie modelu LSTM?

Pracuję z LSTM w Keras / Python, podążając za wieloma przykładami z machinelearningmastery.com , z których oparty jest mój przykładowy kod (poniżej). Zastanawiam się nad przekształceniem problemu w klasyfikację w dyskretne przedziały, ponieważ daje to pewność dla poszczególnych klas, ale wydaje się to złym rozwiązaniem.

Istnieje kilka podobnych tematów (takich jak poniżej), ale wydaje się, że nic nie dotyczy bezpośrednio interwałów prognozowania z sieci neuronowych LSTM (lub rzeczywiście innych):

/stats/25055/how-to-calculate-the-confidence-interval-for-time-series-prediction

Prognozowanie szeregów czasowych przy użyciu ARIMA vs 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)
4Oh4
źródło

Odpowiedzi:

10

Bezpośrednio nie jest to możliwe. Jeśli jednak modelujesz to w inny sposób, możesz uzyskać przedziały ufności. Możesz zamiast normalnej regresji podejść do tego jako oszacowanie ciągłego rozkładu prawdopodobieństwa. Robiąc to na każdym kroku, możesz zaplanować swoją dystrybucję. Sposobami na to są Kernel Mixture Networks ( https://janvdvegt.github.io/2017/06/07/Kernel-Mixture-Networks.html , ujawnienie, mój blog) lub Density Mixture Networks ( http: //www.cedar .buffalo.edu / ~ srihari / CSE574 / Chap5 / Chap5.7-MixDensityNetworks.pdf ), pierwszy wykorzystuje jądra jako bazę i ocenia mieszaninę na tych jądrach, a drugi ocenia mieszaninę rozkładów, w tym parametry każdego z dystrybucje. Prawdopodobieństwa dziennika używa się do szkolenia modelu.

Inną opcją modelowania niepewności jest użycie przerwania podczas treningu, a następnie także podczas wnioskowania. Robisz to wiele razy i za każdym razem, gdy pobierasz próbkę z tylnej części ciała. Nie dostajesz dystrybucji, tylko próbki, ale jest najłatwiejszy do wdrożenia i działa bardzo dobrze.

W twoim przypadku musisz pomyśleć o sposobie generowania t + 2 do t + 10. W zależności od bieżącej konfiguracji może być konieczne próbkowanie z poprzedniego kroku czasowego i karmienie go na następny. To nie działa dobrze z pierwszym podejściem, ani z drugim. Jeśli masz 10 wyjść na krok czasu (t + 1 do t + 10), wszystkie te podejścia są bardziej czyste, ale nieco mniej intuicyjne.

Jan van der Vegt
źródło
2
Ciekawe jest używanie sieci mikstur, spróbuję to zaimplementować. Istnieją pewne stałe badania nad użyciem przerywania tutaj: arxiv.org/abs/1709.01907 i arxiv.org/abs/1506.02142
4Oh4
Uwaga na temat porzucenia, możesz faktycznie obliczyć wariancję przewidywania porzucenia Monte Carlo i użyć tego jako kwantyfikacji niepewności
Charles Chow
To prawda @CharlesChow, ale jest to zły sposób na skonstruowanie przedziału ufności w tym kontekście. Lepiej byłoby posortować wartości i użyć kwantyli ze względu na potencjalnie bardzo wypaczony rozkład.
Jan van der Vegt,
Zgadzam się @JanvanderVegt, ale nadal możesz oszacować statystyki rezygnacji MC bez założenia dystrybucji wyjściowej, to znaczy, że możesz również użyć percentyla lub bootstrapowania do skonstruowania CI rezygnacji MC
Charles Chow
2

Conformal Prediction jako buzz word może być dla Ciebie interesujący, ponieważ działa w wielu warunkach - w szczególności nie wymaga normalnego rozproszonego błędu i działa w prawie każdym modelu uczenia maszynowego.

Dwa miłe wprowadzenia przedstawili Scott Locklin i Henrik Linusson .

Boris W.
źródło
1

Będę się nieco rozbierać i argumentować, że przedział ufności obliczeń jest w praktyce zwykle nieoceniony. Powodem jest to, że zawsze musisz przyjąć cały szereg założeń. Musisz mieć nawet najprostszą regresję liniową

  • Zależność liniowa.
  • Wielowymiarowa normalność.
  • Brak lub niewielka wielokoliniowość.
  • Brak autokorelacji.
  • Homoscedastyczność.

O wiele bardziej pragmatycznym podejściem jest przeprowadzenie symulacji Monte Carlo. Jeśli już wiesz lub chcesz założyć wokół rozkładu zmiennych wejściowych, weź całą próbkę i podaj ją LSTM, teraz możesz empirycznie obliczyć „przedział ufności”.

Louis T.
źródło
1

Tak, możesz. Jedyne, co musisz zmienić, to funkcja utraty. Zaimplementuj funkcję straty stosowaną w regresji kwantowej i zintegruj ją. Ponadto chcesz spojrzeć na to, jak oceniasz te interwały. W tym celu użyłbym wskaźników ICP, MIL i RMIL.

Inigo
źródło