使用LSTM深度學習模型進行溫度的時間序列單步和多步預測

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

本文展示了如何進行多步預測並在模型中使用多個特徵。

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

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

所使用的特徵是過去每小時的溫度數據、每日及每年的循環信號、氣壓及風速。

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

這裡和整篇文章的主數據對象被稱為d。它是通過讀取原始數據創建的:

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)

數據集中共有271008個數據點。

數據似乎是具有明確的周期模式。

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

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

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

我們想讓機器知道，23點和0點比小時0點和4點更接近。我們知道周期是24小時。我們可以用cos(x)和sin(x)函數。函數中的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如下:

新創建的特徵捕捉了周期性模式。 可能會出現一個問題，為什麼我們同時使用sin和cos函數？

在上圖中繪製一條水平線並僅分析其中一條曲線，我們將得到例如cos（7.5h）= cos（17.5h）等。 在學習和預測時，這可能會導致一些錯誤，因此為了使每個點都唯一，我們添加了另一個循環函數。 同時使用這兩個功能，可以將所有時間區分開。

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

# 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"]]

在本節中，我們從datetime列中創建了4個其他功能：daysin，daycos，monthsin和monthcos。

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

要查看溫度與兩個特徵之間的任何關係，我們可以繪製二維直方圖：

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

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

我們使用所有要素工程獲得的數據是：

我們要近似的函數f為：

目標是使用過去的值來預測未來。 數據是時間序列或序列。 對於序列建模，我們將選擇具有LSTM層的遞歸神經網絡的Tensorflow實現。

LSTM網絡的輸入是3D張量：

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

樣本—用於訓練的序列總數。

timesteps-樣本的長度。

功能-使用的功能數量。

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

例如，如果我們假設整個數據是數據的前10行，那麼我們將過去3個小時用作特徵，並希望預測出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

例如，如果我們假設整個數據是數據的前10行，那麼我們將過去3個小時用作特徵，並希望預測出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行數據和7列。 之所以有6個觀測值，是因為前3個滯後被丟棄並且僅用作X數據，並且我們預測提前1步，因此最後一個觀測值也會丟失。

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

最終模型的超參數列表：

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

創建用於建模之前的最後一步是縮放數據。

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

現在我們將數據分為訓練和驗證

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

數據的最終形狀：

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)

剩下的就是使用模型類創建對象，訓練模型並檢查驗證集中的結果。

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

使用訓練好的模型，我們可以預測值並將其與原始值進行比較。

# 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』]]

結果可視化

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

使用訓練好的模型，我們可以預測值並將其與原始值進行比較。

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

要預測提前24小時，唯一需要做的就是更改超參數。 具體來說，是n_ahead變量。 該模型將嘗試使用之前（一周）的168小時來預測接下來的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()

檢查一些隨機的24小時區間：

有些24小時序列似乎彼此接近，而其他序列則不然。

平均絕對誤差為1.69 C，中位數為1.27C。

總結，本文介紹了在對時間序列數據進行建模和預測時使用的簡單管道示例：

讀取，清理和擴充輸入數據

為滯後和n步選擇超參數

為深度學習模型選擇超參數

初始化NNMultistepModel（）類

擬合模型

預測未來n_steps

最後本文的完整代碼：https://github.com/Eligijus112/Vilnius-weather-LSTM

作者：Eligijus Bujokas

deephub翻譯組