【RNN入門到實戰】LSTM從入門到實戰——實現空氣質量預測

語言: CN / TW / HK

摘要

LSTM是一種時間遞迴神經網路,它出現的原因是為了解決RNN的一個致命的缺陷。RNN在處理長期依賴(時間序列上距離較遠的節點)時,因為計算距離較遠的節點之間的聯絡時會涉及雅可比矩陣的多次相乘,會造成梯度消失或者梯度膨脹的現象。為了解決該問題,研究人員提出了許多解決辦法,例如ESN(Echo State Network),增加有漏單元(Leaky Units)等等。其中最成功應用最廣泛的就是門限RNN(Gated RNN),而LSTM就是門限RNN中最著名的一種。有漏單元通過設計連線間的權重係數,從而允許RNN累積距離較遠節點間的長期聯絡;而門限RNN則泛化了這樣的思想,允許在不同時刻改變該係數,且允許網路忘記當前已經累積的資訊。\ 在這裡插入圖片描述

RNN和LSTM的區別

所有 RNN 都具有一種重複神經網路模組的鏈式的形式。在標準的 RNN 中,這個重複的模組只有一個非常簡單的結構,例如一個 tanh 層,如下圖所示:\ 在這裡插入圖片描述LSTM 同樣是這樣的結構,但是重複的模組擁有一個不同的結構。不同於單一神經網路層,這裡是有四個,以一種非常特殊的方式進行互動。\ LSTM結構

詳解LSTM

LSTM 的關鍵就是Cell狀態,水平線在圖上方貫穿執行。Cell狀態類似於傳送帶。直接在整個鏈上執行,只有一些少量的線性互動。資訊在上面流傳保持不變會很容易。示意圖如下所示:\ 在這裡插入圖片描述

LSTM 有通過精心設計的稱作為“門”的結構來去除或者增加資訊到Cell狀態的能力。門是一種讓資訊選擇式通過的方法。他們包含一個 sigmoid 神經網路層和一個 pointwise 乘法操作。示意圖如下:

在這裡插入圖片描述

\ LSTM 擁有三個門,分別是忘記層門,輸入層門和輸出層門,來保護和控制Cell狀態。\ 忘記層門\ ​ 作用物件:Cell狀態 。\ ​ 作用:將細胞狀態中的資訊選擇性的遺忘。\ ​ 操作步驟:該門會讀取 h t − 1 h_{t-1} ht−1​和 x t x_t xt​,輸出一個在 0 到 1 之間的數值給每個在細胞狀態 C t − 1 C_{t-1} Ct−1​中的數字。1 表示“完全保留”,0 表示“完全捨棄”。示意圖如下:\ 在這裡插入圖片描述

輸入層門\ ​ 作用物件:細胞狀態\ ​ 作用:將新的資訊選擇性的記錄到細胞狀態中。\ ​ 操作步驟:\ ​ 步驟一,sigmoid 層稱 “輸入門層” 決定什麼值我們將要更新。\ ​ 步驟二,tanh 層建立一個新的候選值向量 C \~ t \tilde{C}_t C\~t​加入到狀態中。其示意圖如下:\ 在這裡插入圖片描述\ ​ 步驟三:將 c t − 1 c_{t-1} ct−1​更新為 c t c_{t} ct​。將舊狀態與 f t f_t ft​相乘,丟棄掉我們確定需要丟棄的資訊。接著加上 i t ∗ C \~ t i_t * \tilde{C}_t it​∗C\~t​得到新的候選值,根據我們決定更新每個狀態的程度進行變化。其示意圖如下:

在這裡插入圖片描述\ 動圖演示\ 在這裡插入圖片描述

輸出層門\ ​ 作用物件:隱層 h t h_t ht​\ ​ 作用:確定輸出什麼值。\ ​ 操作步驟:\ ​ 步驟一:通過sigmoid 層來確定細胞狀態的哪個部分將輸出。\ ​ 步驟二:把細胞狀態通過 tanh 進行處理,並將它和 sigmoid 門的輸出相乘,最終我們僅僅會輸出我們確定輸出的那部分。\ 其示意圖如下所示:\ 在這裡插入圖片描述

動圖演示\ 在這裡插入圖片描述

實戰——使用LSTM實現空氣質量預測

資料來源自位於北京的美國大使館在2010年至2014年共5年間每小時採集的天氣及空氣汙染指數。\   資料集包括日期、PM2.5濃度、露點、溫度、風向、風速、累積小時雪量和累積小時雨量。原始資料中完整的特徵如下:\ 1.No 行數\ 2.year 年\ 3.month 月\ 4.day 日\ 5.hour 小時\ 6.pm2.5 PM2.5濃度\ 7.DEWP 露點\ 8.TEMP 溫度\ 9.PRES 大氣壓\ 10.cbwd 風向\ 11.lws 風速\ 12.ls 累積雪量\ 13.lr 累積雨量\ 我們可以利用此資料集搭建預測模型,利用前一個或幾個小時的天氣條件和汙染資料預測下一個(當前)時刻的汙染程度。

資料處理

首先,我們必須清洗資料。以下是原始資料集的前幾行。

python No year month day hour pm2.5 DEWP TEMP PRES cbwd Iws Is Ir 0 1 2010 1 1 0 NaN -21 -11.0 1021.0 NW 1.79 0 0 1 2 2010 1 1 1 NaN -21 -12.0 1020.0 NW 4.92 0 0 2 3 2010 1 1 2 NaN -21 -11.0 1019.0 NW 6.71 0 0 3 4 2010 1 1 3 NaN -21 -14.0 1019.0 NW 9.84 0 0 4 5 2010 1 1 4 NaN -20 -12.0 1018.0 NW 12.97 0 0 5 6 2010 1 1 5 NaN -19 -10.0 1017.0 NW 16.10 0 0 6 7 2010 1 1 6 NaN -19 -9.0 1017.0 NW 19.23 0 0 7 8 2010 1 1 7 NaN -19 -9.0 1017.0 NW 21.02 0 0 8 9 2010 1 1 8 NaN -19 -9.0 1017.0 NW 24.15 0 0 9 10 2010 1 1 9 NaN -20 -8.0 1017.0 NW 27.28 0 0

資料理清的步驟:\   1、將year, month, day, hour四列整合為一個日期時間。\   2、刪除No列,這個列對於資料預測沒有作用,如果有作用說明見鬼了。\   3、將資料集中所有的NaN設定為0,NaN沒有辦法用來計算。\   4、刪除前24行,前24行的pm2.5沒有記錄,留著沒有用。\ 完整的程式碼如下:

```python

from pandas import read_csv from datetime import datetime

load data

def parse(x): return datetime.strptime(x, '%Y %m %d %H')

讀取資料,將year, month, day, hour四列合併成一列。

dataset = read_csv('raw.csv', parse_dates = [['year', 'month', 'day', 'hour']], index_col=0, date_parser=parse)

刪除No列

dataset.drop('No', axis=1, inplace=True)

修改列名

dataset.columns = ['pollution', 'dew', 'temp', 'press', 'wnd_dir', 'wnd_spd', 'snow', 'rain'] dataset.index.name = 'date' print(dataset)

將所有的NaN設定為0

dataset['pollution'].fillna(0, inplace=True)

刪除前24行

dataset = dataset[24:]

瀏覽前5行資料

print(dataset.head(5))

save to file

dataset.to_csv('pollution.csv') ```

載入了“pollution.csv”檔案,並對除了類別型特性“風速”的每一列資料分別繪圖。

python dataset = pd.read_csv('pollution.csv', header=0, index_col=0) values = dataset.values # specify columns to plot groups = [0, 1, 2, 3, 5, 6, 7] i = 1 # plot each column pyplot.figure(figsize=(10, 10)) for group in groups: pyplot.subplot(len(groups), 1, i) pyplot.plot(values[:, group]) pyplot.title(dataset.columns[group], y=0.5, loc='right') i += 1 pyplot.show()

執行上面的程式碼,並對7個變數在5年的範圍內繪圖。 在這裡插入圖片描述\ 利用sklearn的預處理模組對類別特徵“風向”進行編碼,當然也可以對該特徵進行one-hot編碼。 接著對所有的特徵進行歸一化處理,然後將資料集轉化為有監督學習問題,同時將需要預測的當前時刻(t)的天氣條件特徵移除,程式碼如下:

```python def series_to_supervised(data, n_in=1, n_out=1, dropnan=True): # convert series to supervised learning n_vars = 1 if type(data) is list else data.shape[1] df = pd.DataFrame(data) cols, names = list(), list() # input sequence (t-n, ... t-1) for i in range(n_in, 0, -1): cols.append(df.shift(i)) names += [('var%d(t-%d)' % (j + 1, i)) for j in range(n_vars)] # forecast sequence (t, t+1, ... t+n) for i in range(0, n_out): cols.append(df.shift(-i)) if i == 0: names += [('var%d(t)' % (j + 1)) for j in range(n_vars)] else: names += [('var%d(t+%d)' % (j + 1, i)) for j in range(n_vars)] # put it all together agg = pd.concat(cols, axis=1) agg.columns = names # drop rows with NaN values if dropnan: agg.dropna(inplace=True) return agg

load dataset

dataset = pd.read_csv('pollution.csv', header=0, index_col=0) values = dataset.values

integer encode direction

encoder = LabelEncoder() print(values[:, 4]) values[:, 4] = encoder.fit_transform(values[:, 4]) print(values[:, 4])

ensure all data is float

values = values.astype('float32')

normalize features

scaler = MinMaxScaler(feature_range=(0, 1)) scaled = scaler.fit_transform(values)

frame as supervised learning

reframed = series_to_supervised(scaled, 1, 1)

drop columns we don't want to predict

reframed.drop(reframed.columns[[9, 10, 11, 12, 13, 14, 15]], axis=1, inplace=True) print(reframed.head()) ```

構造模型

首先,我們需要將處理後的資料集劃分為訓練集和測試集。為了加速模型的訓練,我們僅利用第一年資料進行訓練,然後利用剩下的4年進行評估。\   下面的程式碼將資料集進行劃分,然後將訓練集和測試集劃分為輸入和輸出變數,最終將輸入(X)改造為LSTM的輸入格式,即[samples,timesteps,features]。

```python

split into train and test sets

values = reframed.values n_train_hours = 365 * 24 train = values[:n_train_hours, :] test = values[n_train_hours:, :]

split into input and outputs

train_X, train_y = train[:, :-1], train[:, -1] test_X, test_y = test[:, :-1], test[:, -1]

reshape input to be 3D [samples, timesteps, features]

train_X = train_X.reshape((train_X.shape[0], 1, train_X.shape[1])) test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1])) print(train_X.shape, train_y.shape, test_X.shape, test_y.shape) ```

執行上述程式碼列印訓練集和測試集的輸入輸出格式,其中9K小時資料作訓練集,35K小時資料作測試集。\ (8760, 1, 8) (8760,) (35039, 1, 8) (35039,)\ 現在可以搭建LSTM模型了。 LSTM模型中,隱藏層有50個神經元,輸出層1個神經元(迴歸問題),輸入變數是一個時間步(t-1)的特徵,損失函式採用Mean Absolute Error(MAE),優化演算法採用Adam,模型採用50個epochs並且每個batch的大小為72。\   最後,在fit()函式中設定validation_data引數,記錄訓練集和測試集的損失,並在完成訓練和測試後繪製損失圖。

```python checkpointer = ModelCheckpoint(filepath='best_model.hdf5', monitor='val_loss', verbose=1, save_best_only=True, mode='min') reduce = ReduceLROnPlateau(monitor='val_loss', patience=10, verbose=1, factor=0.5, min_lr=1e-6)

model = Sequential() model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2]))) model.add(Dense(1)) model.compile(loss='mae', optimizer='adam') # fit network history = model.fit(train_X, train_y, epochs=300, batch_size=64, validation_data=(test_X, test_y), verbose=1, callbacks=[checkpointer, reduce], shuffle=True) # plot history pyplot.plot(history.history['loss'], label='train') pyplot.plot(history.history['val_loss'], label='test') pyplot.legend() pyplot.show() ```

模型評估

接下里我們對模型效果進行評估。\   值得注意的是:需要將預測結果和部分測試集資料組合然後進行比例反轉(invert the scaling),同時也需要將測試集上的預期值也進行比例轉換。\   (We combine the forecast with the test dataset and invert the scaling. We also invert scaling on the test dataset with the expected pollution numbers.)\   至於在這裡為什麼進行比例反轉,是因為我們將原始資料進行了預處理(連同輸出值y),此時的誤差損失計算是在處理之後的資料上進行的,為了計算在原始比例上的誤差需要將資料進行轉化。同時筆者有個小Tips:就是反轉時的矩陣大小一定要和原來的大小(shape)完全相同,否則就會報錯。\   通過以上處理之後,再結合RMSE(均方根誤差)計算損失。

```python yhat = model.predict(test_X) test_X = test_X.reshape((test_X.shape[0], test_X.shape[2]))

invert scaling for forecast

inv_yhat = concatenate((yhat, test_X[:, 1:]), axis=1) inv_yhat = scaler.inverse_transform(inv_yhat) inv_yhat = inv_yhat[:, 0]

invert scaling for actual

inv_y = scaler.inverse_transform(test_X) inv_y = inv_y[:, 0]

calculate RMSE

rmse = sqrt(mean_squared_error(inv_y, inv_yhat)) print('Test RMSE: %.3f' % rmse) ```

完整程式碼

```python import pandas as pd from datetime import datetime from matplotlib import pyplot from sklearn.preprocessing import LabelEncoder, MinMaxScaler from sklearn.metrics import mean_squared_error from tensorflow.keras.models import Sequential from tensorflow.keras.layers import Dense from tensorflow.keras.layers import LSTM from numpy import concatenate from math import sqrt

load data

def parse(x): return datetime.strptime(x, '%Y %m %d %H')

def read_raw(): dataset = pd.read_csv('raw.csv', parse_dates=[['year', 'month', 'day', 'hour']], index_col=0, date_parser=parse) dataset.drop('No', axis=1, inplace=True) # manually specify column names dataset.columns = ['pollution', 'dew', 'temp', 'press', 'wnd_dir', 'wnd_spd', 'snow', 'rain'] dataset.index.name = 'date' # mark all NA values with 0 dataset['pollution'].fillna(0, inplace=True) # drop the first 24 hours dataset = dataset[24:] # summarize first 5 rows print(dataset.head(5)) # save to file dataset.to_csv('pollution.csv')

def drow_pollution(): dataset = pd.read_csv('pollution.csv', header=0, index_col=0) values = dataset.values # specify columns to plot groups = [0, 1, 2, 3, 5, 6, 7] i = 1 # plot each column pyplot.figure(figsize=(10, 10)) for group in groups: pyplot.subplot(len(groups), 1, i) pyplot.plot(values[:, group]) pyplot.title(dataset.columns[group], y=0.5, loc='right') i += 1 pyplot.show()

def series_to_supervised(data, n_in=1, n_out=1, dropnan=True): # convert series to supervised learning n_vars = 1 if type(data) is list else data.shape[1] df = pd.DataFrame(data) cols, names = list(), list() # input sequence (t-n, ... t-1) for i in range(n_in, 0, -1): cols.append(df.shift(i)) names += [('var%d(t-%d)' % (j + 1, i)) for j in range(n_vars)] # forecast sequence (t, t+1, ... t+n) for i in range(0, n_out): cols.append(df.shift(-i)) if i == 0: names += [('var%d(t)' % (j + 1)) for j in range(n_vars)] else: names += [('var%d(t+%d)' % (j + 1, i)) for j in range(n_vars)] # put it all together agg = pd.concat(cols, axis=1) agg.columns = names # drop rows with NaN values if dropnan: agg.dropna(inplace=True) return agg

def cs_to_sl(): # load dataset dataset = pd.read_csv('pollution.csv', header=0, index_col=0) values = dataset.values # integer encode direction encoder = LabelEncoder() print(values[:, 4]) values[:, 4] = encoder.fit_transform(values[:, 4]) print(values[:, 4]) # ensure all data is float values = values.astype('float32') # normalize features scaler = MinMaxScaler(feature_range=(0, 1)) scaled = scaler.fit_transform(values) # frame as supervised learning reframed = series_to_supervised(scaled, 1, 1) # drop columns we don't want to predict reframed.drop(reframed.columns[[9, 10, 11, 12, 13, 14, 15]], axis=1, inplace=True) print(reframed.head()) return reframed, scaler

def train_test(reframed): # split into train and test sets values = reframed.values n_train_hours = 365 * 24 train = values[:n_train_hours, :] test = values[n_train_hours:, :] # split into input and outputs train_X, train_y = train[:, :-1], train[:, -1] test_X, test_y = test[:, :-1], test[:, -1] # reshape input to be 3D [samples, timesteps, features] train_X = train_X.reshape((train_X.shape[0], 1, train_X.shape[1])) test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1])) print(train_X.shape, train_y.shape, test_X.shape, test_y.shape) return train_X, train_y, test_X, test_y

def fit_network(train_X, train_y, test_X, test_y, scaler): print(train_X.shape) print(train_X.shape[1]) model = Sequential() model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2]))) model.add(Dense(1)) model.compile(loss='mae', optimizer='adam') # fit network history = model.fit(train_X, train_y, epochs=50, batch_size=72, validation_data=(test_X, test_y), verbose=2, shuffle=False) # plot history pyplot.plot(history.history['loss'], label='train') pyplot.plot(history.history['val_loss'], label='test') pyplot.legend() pyplot.show() # make a prediction yhat = model.predict(test_X) test_X = test_X.reshape((test_X.shape[0], test_X.shape[2])) # invert scaling for forecast inv_yhat = concatenate((yhat, test_X[:, 1:]), axis=1) inv_yhat = scaler.inverse_transform(inv_yhat) inv_yhat = inv_yhat[:, 0] # invert scaling for actual inv_y = scaler.inverse_transform(test_X) inv_y = inv_y[:, 0] # calculate RMSE rmse = sqrt(mean_squared_error(inv_y, inv_yhat)) print('Test RMSE: %.3f' % rmse)

if name == 'main': drow_pollution() reframed, scaler = cs_to_sl() train_X, train_y, test_X, test_y = train_test(reframed) fit_network(train_X, train_y, test_X, test_y, scaler) ```

程式碼和資料集連結:https\://download.csdn.net/download/hhhhhhhhhhwwwwwwwwww/19781047