使用LSTM深度學(xué)習(xí)模型進(jìn)行溫度的時間序列單步和多步預(yù)測

禁忌石 2022-04-10

展開全文

本文的目的是提供代碼示例，并解釋使用python和TensorFlow建模時間序列數(shù)據(jù)的思路。

本文展示了如何進(jìn)行多步預(yù)測并在模型中使用多個特征。

本文的簡單版本是，使用過去48小時的數(shù)據(jù)和對未來1小時的預(yù)測(一步)，我獲得了溫度誤差的平均絕對誤差0.48(中值0.34)度。

利用過去168小時的數(shù)據(jù)并提前24小時進(jìn)行預(yù)測，平均絕對誤差為攝氏溫度1.69度(中值1.27)。

所使用的特征是過去每小時的溫度數(shù)據(jù)、每日及每年的循環(huán)信號、氣壓及風(fēng)速。

使用來自
https:///的API獲取數(shù)據(jù)。這些數(shù)據(jù)從1990年1月1日到2020.11月30日每小時在維爾紐斯電視塔附近收集一次。維爾紐斯不是一個大城市，電視塔就在城市里，所以電視塔附近的溫度應(yīng)該和城市所有地方的溫度非常相似。

這里和整篇文章的主數(shù)據(jù)對象被稱為d。它是通過讀取原始數(shù)據(jù)創(chuàng)建的:

d = pd.read_csv('data/weather.csv’)# Converting the dt column to datetime object d['dt’] = [datetime.datetime.utcfromtimestamp(x) for x in d['dt’]]# Sorting by the date d.sort_values('dt’, inplace=True)

數(shù)據(jù)集中共有271008個數(shù)據(jù)點。

數(shù)據(jù)似乎是具有明確的周期模式。

上面的圖表顯示，氣溫有一個清晰的晝夜循環(huán)——中間溫度在中午左右最高，在午夜左右最低。

這種循環(huán)模式在按月份分組的溫度上更為明顯——最熱的月份是6月到8月，最冷的月份是12月到2月。

數(shù)據(jù)現(xiàn)在的問題是，我們只有date列。如果將其轉(zhuǎn)換為數(shù)值(例如，提取時間戳(以秒為單位))并將其作為建模時的特性添加，那么循環(huán)特性將丟失。因此，我們需要做的第一件事就是設(shè)計一些能夠抓住周期性趨勢的特性。

我們想讓機器知道，23點和0點比小時0點和4點更接近。我們知道周期是24小時。我們可以用cos(x)和sin(x)函數(shù)。函數(shù)中的x是一天中的一個小時。

# Extracting the hour of dayd['hour'] = [x.hour for x in d['dt']]# Creating the cyclical daily feature d['day_cos'] = [np.cos(x * (2 * np.pi / 24)) for x in d['hour']]d['day_sin'] = [np.sin(x * (2 * np.pi / 24)) for x in d['hour']]

得到的dataframe如下:

新創(chuàng)建的特征捕捉了周期性模式。可能會出現(xiàn)一個問題，為什么我們同時使用sin和cos函數(shù)？

在上圖中繪制一條水平線并僅分析其中一條曲線，我們將得到例如cos（7.5h）= cos（17.5h）等。在學(xué)習(xí)和預(yù)測時，這可能會導(dǎo)致一些錯誤，因此為了使每個點都唯一，我們添加了另一個循環(huán)函數(shù)。同時使用這兩個功能，可以將所有時間區(qū)分開。

為了在一年中的某個時間創(chuàng)建相同的循環(huán)邏輯，我們將使用時間戳功能。 python中的時間戳是一個值，用于計算自1970.01.01 0H：0m：0s以來經(jīng)過了多少秒。 python中的每個date對象都具有timestamp（）函數(shù)。

# Extracting the timestamp from the datetime object d['timestamp'] = [x.timestamp() for x in d['dt']]# Seconds in day s = 24 * 60 * 60# Seconds in year year = (365.25) * sd['month_cos'] = [np.cos((x) * (2 * np.pi / year)) for x in d['timestamp']]d['month_sin'] = [np.sin((x) * (2 * np.pi / year)) for x in d['timestamp']]

在本節(jié)中，我們從datetime列中創(chuàng)建了4個其他功能：daysin，daycos，monthsin和monthcos。

在天氣數(shù)據(jù)集中，還有兩列：wind_speed和pressure。風(fēng)速以米/秒（m / s）為單位，壓力以百帕斯卡（hPa）為單位。

要查看溫度與兩個特征之間的任何關(guān)系，我們可以繪制二維直方圖：

顏色越強烈，兩個分布的某些bin值之間的關(guān)系就越大。例如，當(dāng)壓力在1010和1020 hPa左右時，溫度往往會更高。

我們還將在建模中使用這兩個功能。

我們使用所有要素工程獲得的數(shù)據(jù)是：

我們要近似的函數(shù)f為：

目標(biāo)是使用過去的值來預(yù)測未來。數(shù)據(jù)是時間序列或序列。對于序列建模，我們將選擇具有LSTM層的遞歸神經(jīng)網(wǎng)絡(luò)的Tensorflow實現(xiàn)。

LSTM網(wǎng)絡(luò)的輸入是3D張量：

（樣本，時間步長，功能）

樣本—用于訓(xùn)練的序列總數(shù)。

timesteps-樣本的長度。

功能-使用的功能數(shù)量。

建模之前的第一件事是將2D格式的數(shù)據(jù)轉(zhuǎn)換為3D數(shù)組。以下功能可以做到這一點：

例如，如果我們假設(shè)整個數(shù)據(jù)是數(shù)據(jù)的前10行，那么我們將過去3個小時用作特征，并希望預(yù)測出1步：

def create_X_Y(ts: np.array, lag=1, n_ahead=1, target_index=0) -> tuple:    '''    A method to create X and Y matrix from a time series array for the training of     deep learning models     '''    # Extracting the number of features that are passed from the array     n_features = ts.shape[1]        # Creating placeholder lists    X, Y = [], []    if len(ts) - lag <= 0:        X.append(ts)    else:        for i in range(len(ts) - lag - n_ahead):            Y.append(ts[(i + lag):(i + lag + n_ahead), target_index])            X.append(ts[i:(i + lag)])    X, Y = np.array(X), np.array(Y)    # Reshaping the X array to an RNN input shape     X = np.reshape(X, (X.shape[0], lag, n_features))    return X, Y

例如，如果我們假設(shè)整個數(shù)據(jù)是數(shù)據(jù)的前10行，那么我們將過去3個小時用作特征，并希望預(yù)測出1步：

ts = d['temp’, 'day_cos’, 'day_sin’, 'month_sin’, 'month_cos’, 'pressure’, 'wind_speed’].head(10).valuesX, Y = create_X_Y(ts, lag=3, n_ahead=1)

如我們所見，X矩陣的形狀是6個樣本，3個時間步長和7個特征。換句話說，我們有6個觀測值，每個觀測值都有3行數(shù)據(jù)和7列。之所以有6個觀測值，是因為前3個滯后被丟棄并且僅用作X數(shù)據(jù)，并且我們預(yù)測提前1步，因此最后一個觀測值也會丟失。

上圖中顯示了X和Y的第一個值對。

最終模型的超參數(shù)列表：

# Number of lags (hours back) to use for modelslag = 48# Steps ahead to forecast n_ahead = 1# Share of obs in testing test_share = 0.1# Epochs for trainingepochs = 20# Batch size batch_size = 512# Learning ratelr = 0.001# Number of neurons in LSTM layern_layer = 10# The features used in the modeling features_final = ['temp’, 'day_cos’, 'day_sin’, 'month_sin’, 'month_cos’, 'pressure’, 'wind_speed’]

模型代碼

class NNMultistepModel():        def __init__(        self,         X,         Y,         n_outputs,        n_lag,        n_ft,        n_layer,        batch,        epochs,         lr,        Xval=None,        Yval=None,        mask_value=-999.0,        min_delta=0.001,        patience=5    ):        lstm_input = Input(shape=(n_lag, n_ft))        # Series signal         lstm_layer = LSTM(n_layer, activation='relu')(lstm_input)        x = Dense(n_outputs)(lstm_layer)                self.model = Model(inputs=lstm_input, outputs=x)        self.batch = batch         self.epochs = epochs        self.n_layer=n_layer        self.lr = lr         self.Xval = Xval        self.Yval = Yval        self.X = X        self.Y = Y        self.mask_value = mask_value        self.min_delta = min_delta        self.patience = patience    def trainCallback(self):        return EarlyStopping(monitor='loss', patience=self.patience, min_delta=self.min_delta)    def train(self):        # Getting the untrained model         empty_model = self.model                # Initiating the optimizer        optimizer = keras.optimizers.Adam(learning_rate=self.lr)        # Compiling the model        empty_model.compile(loss=losses.MeanAbsoluteError(), optimizer=optimizer)        if (self.Xval is not None) & (self.Yval is not None):            history = empty_model.fit(                self.X,                 self.Y,                 epochs=self.epochs,                 batch_size=self.batch,                 validation_data=(self.Xval, self.Yval),                 shuffle=False,                callbacks=[self.trainCallback()]            )        else:            history = empty_model.fit(                self.X,                 self.Y,                 epochs=self.epochs,                 batch_size=self.batch,                shuffle=False,                callbacks=[self.trainCallback()]            )                # Saving to original model attribute in the class        self.model = empty_model                # Returning the training history        return history        def predict(self, X):        return self.model.predict(X)

創(chuàng)建用于建模之前的最后一步是縮放數(shù)據(jù)。

# Subseting only the needed columns ts = d[features_final]nrows = ts.shape[0]# Spliting into train and test setstrain = ts[0:int(nrows * (1 — test_share))]test = ts[int(nrows * (1 — test_share)):]# Scaling the data train_mean = train.mean()train_std = train.std()train = (train — train_mean) / train_stdtest = (test — train_mean) / train_std# Creating the final scaled frame ts_s = pd.concat([train, test])# Creating the X and Y for trainingX, Y = create_X_Y(ts_s.values, lag=lag, n_ahead=n_ahead)n_ft = X.shape[2]

現(xiàn)在我們將數(shù)據(jù)分為訓(xùn)練和驗證

# Spliting into train and test sets Xtrain, Ytrain = X[0:int(X.shape[0] * (1 — test_share))], Y[0:int(X.shape[0] * (1 — test_share))]Xval, Yval = X[int(X.shape[0] * (1 — test_share)):], Y[int(X.shape[0] * (1 — test_share)):]

數(shù)據(jù)的最終形狀：

Shape of training data: (243863, 48, 7)Shape of the target data: (243863, 1)Shape of validation data: (27096, 48, 7)Shape of the validation target data: (27096, 1)

剩下的就是使用模型類創(chuàng)建對象，訓(xùn)練模型并檢查驗證集中的結(jié)果。

# Initiating the model objectmodel = NNMultistepModel( X=Xtrain, Y=Ytrain, n_outputs=n_ahead, n_lag=lag, n_ft=n_ft, n_layer=n_layer, batch=batch_size, epochs=epochs,  lr=lr, Xval=Xval, Yval=Yval,)# Training of the model history = model.train()

使用訓(xùn)練好的模型，我們可以預(yù)測值并將其與原始值進(jìn)行比較。

# Comparing the forecasts with the actual valuesyhat = [x[0] for x in model.predict(Xval)]y = [y[0] for y in Yval]# Creating the frame to store both predictionsdays = d['dt’].values[-len(y):]frame = pd.concat([ pd.DataFrame({'day’: days, 'temp’: y, 'type’: 'original’}), pd.DataFrame({'day’: days, 'temp’: yhat, 'type’: 'forecast’})])# Creating the unscaled values columnframe['temp_absolute’] = [(x * train_std['temp’]) + train_mean['temp’] for x in frame['temp’]]# Pivotingpivoted = frame.pivot_table(index=’day’, columns=’type’)pivoted.columns = ['_’.join(x).strip() for x in pivoted.columns.values]pivoted['res’] = pivoted['temp_absolute_original’] — pivoted['temp_absolute_forecast’]pivoted['res_abs’] = [abs(x) for x in pivoted['res’]]

結(jié)果可視化

plt.figure(figsize=(12, 12))plt.plot(pivoted.index, pivoted.temp_absolute_original, color=’blue’, label=’original’)plt.plot(pivoted.index, pivoted.temp_absolute_forecast, color=’red’, label=’forecast’, alpha=0.6)plt.title('Temperature forecasts — absolute data’)plt.legend()plt.show()

使用訓(xùn)練好的模型，我們可以預(yù)測值并將其與原始值進(jìn)行比較。

中位數(shù)絕對誤差為0.34攝氏度，平均值為0.48攝氏度。

要預(yù)測提前24小時，唯一需要做的就是更改超參數(shù)。具體來說，是n_ahead變量。該模型將嘗試使用之前（一周）的168小時來預(yù)測接下來的24小時值。

# Number of lags (hours back) to use for modelslag = 168# Steps ahead to forecast n_ahead = 24# Share of obs in testing test_share = 0.1# Epochs for trainingepochs = 20# Batch size batch_size = 512# Learning ratelr = 0.001# Number of neurons in LSTM layern_layer = 10# Creating the X and Y for trainingX, Y = create_X_Y(ts_s.values, lag=lag, n_ahead=n_ahead)n_ft = X.shape[2]Xtrain, Ytrain = X[0:int(X.shape[0] * (1 - test_share))], Y[0:int(X.shape[0] * (1 - test_share))]Xval, Yval = X[int(X.shape[0] * (1 - test_share)):], Y[int(X.shape[0] * (1 - test_share)):]# Creating the model object model = NNMultistepModel(    X=Xtrain,    Y=Ytrain,    n_outputs=n_ahead,    n_lag=lag,    n_ft=n_ft,    n_layer=n_layer,    batch=batch_size,    epochs=epochs,     lr=lr,    Xval=Xval,    Yval=Yval,)# Training the model history = model.train()