순환 신경망 예제 - LSTM 모델을 활용한 베이징 공기질 예측
이번 포스팅에서는 지난번 설명한 순환 신경망 모델 중 LSTM모델을 사용하여 베이징 공기 질 예측을 시도해 본 내용에 대해 다루도록 하겠습니다.
순환 신경망에 대한 내용은 이전 포스팅을 참고하면 됩니다.
https://epozen-dt.github.io/2022/01/26/recurrent-neural-network.html
1. 데이터셋
우선 이번 예제에서 사용한 데이터셋에 대해 설명하도록 하겠습니다.
이번 예제에서 사용한 데이터셋은 Beijing Multi-Site Air-Quality Data Dataset으로, UCI에서 제공하는 데이터셋을 사용하여 진행했습니다.
해당 데이터셋은 다음 링크에서 받아 사용할 수 있습니다.
https://archive.ics.uci.edu/ml/datasets/Beijing+Multi-Site+Air-Quality+Data
데이터셋에 대해 간략하게 소개해보자면,
우선 해당 내용은 베이징의 12 위치에서 시간별로 총 18개의 칼럼으로 구성된 데이터셋입니다.
이 중, 인덱스, 연도, 달, 날짜, 시간을 제외한 13개의 컬럼에 대한 내용은 다음과 같습니다.
| 변수명 | 내용 |
|---|---|
| PM2.5 (target) | PM2.5 농도(ug/m^3 ) |
| PM10 | PM10 농도(ug/m^3) |
| SO2 | SO2 농도(ug/m^3) |
| NO2 | NO2 농도(ug/m^3) |
| CO | CO 농도(ug/m^3) |
| O3 | O3 농도( ug/m^3) |
| TEMP | 온도(섭씨) |
| PRES | 기압(hPa) |
| DEWP | 이슬점 온도(섭씨) |
| RAIN | 강수량(mm) |
| wd | 풍향 |
| WSPM | 풍속(m/s) |
| station | 측정 위치 명 |
2. 예제 수행
예제는 keras를 사용했습니다.
1. 데이터 처리
우선 csv 형식으로된 데이터를 모두 불러오겠습니다.
path= './data/PRSA_Data/'
all_files = glob.glob(path + '/*.csv')
df = pd.DataFrame()
_list = []
for _file in all_files:
_df = pd.read_csv(_file, index_col=None, header=0)
_list.append(_df)
df = pd.concat(_list)
불러온 데이터를 인덱스 제외하고, station과 연도를 기준으로 정렬합니다.
cols = ['year', 'month', 'day', 'hour',
'PM2.5', 'PM10', 'SO2', 'NO2', 'CO', 'O3',
'TEMP', 'PRES', 'DEWP', 'RAIN', 'wd', 'WSPM', 'station'] # NO 제외
df = df[cols]
df = df.sort_values(['station','year'])
보기 편하도록 연도, 달, 일, 시간 칼럼을 묶어 date칼럼으로 만듭니다.
df['date'] = pd.to_datetime(df[['year','month','day','hour']])
df.drop(['year','month','day','hour'],axis=1,inplace=True)
df=df.set_index('date')
df
| PM2.5 | PM10 | SO2 | NO2 | CO | O3 | TEMP | PRES | DEWP | RAIN | wd | WSPM | station | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| date | |||||||||||||
| 2013-03-01 00:00:00 | 4.0 | 4.0 | 4.0 | 7.0 | 300.0 | 77.0 | -0.7 | 1023.0 | -18.8 | 0.0 | NNW | 4.4 | Aotizhongxin |
| 2013-03-01 01:00:00 | 8.0 | 8.0 | 4.0 | 7.0 | 300.0 | 77.0 | -1.1 | 1023.2 | -18.2 | 0.0 | N | 4.7 | Aotizhongxin |
| 2013-03-01 02:00:00 | 7.0 | 7.0 | 5.0 | 10.0 | 300.0 | 73.0 | -1.1 | 1023.5 | -18.2 | 0.0 | NNW | 5.6 | Aotizhongxin |
| 2013-03-01 03:00:00 | 6.0 | 6.0 | 11.0 | 11.0 | 300.0 | 72.0 | -1.4 | 1024.5 | -19.4 | 0.0 | NW | 3.1 | Aotizhongxin |
| 2013-03-01 04:00:00 | 3.0 | 3.0 | 12.0 | 12.0 | 300.0 | 72.0 | -2.0 | 1025.2 | -19.5 | 0.0 | N | 2.0 | Aotizhongxin |
| … | … | … | … | … | … | … | … | … | … | … | … | … | … |
| 2017-02-28 19:00:00 | 11.0 | 32.0 | 3.0 | 24.0 | 400.0 | 72.0 | 12.5 | 1013.5 | -16.2 | 0.0 | NW | 2.4 | Wanshouxigong |
| 2017-02-28 20:00:00 | 13.0 | 32.0 | 3.0 | 41.0 | 500.0 | 50.0 | 11.6 | 1013.6 | -15.1 | 0.0 | WNW | 0.9 | Wanshouxigong |
| 2017-02-28 21:00:00 | 14.0 | 28.0 | 4.0 | 38.0 | 500.0 | 54.0 | 10.8 | 1014.2 | -13.3 | 0.0 | NW | 1.1 | Wanshouxigong |
| 2017-02-28 22:00:00 | 12.0 | 23.0 | 4.0 | 30.0 | 400.0 | 59.0 | 10.5 | 1014.4 | -12.9 | 0.0 | NNW | 1.2 | Wanshouxigong |
| 2017-02-28 23:00:00 | 13.0 | 19.0 | 4.0 | 38.0 | 600.0 | 49.0 | 8.6 | 1014.1 | -15.9 | 0.0 | NNE | 1.3 | Wanshouxigong |
420768 rows × 13 columns
2. Null값 처리
그 다음 데이터셋 중간중간 포함되어있는 결측치를 찾아 처리해줍니다.
이를 처리하는 방법에는 여러가지가 있는데, 이번에는 중간 값을 집어넣어 주는 방식을 사용하겠습니다.
결측치(null) 비율 확인
print(round(df.isnull().sum()/len(df.index), 2)*100)
PM2.5 2.0
PM10 2.0
SO2 2.0
NO2 3.0
CO 5.0
O3 3.0
TEMP 0.0
PRES 0.0
DEWP 0.0
RAIN 0.0
wd 0.0
WSPM 0.0
station 0.0
dtype: float64
각 항목에 대해 중간값으로 처리
df['PM2.5'].fillna(df['PM2.5'].median(), inplace=True)
df['PM10'].fillna(df['PM10'].median(), inplace=True)
df['SO2'].fillna(df['SO2'].median(), inplace=True)
df['NO2'].fillna(df['NO2'].mean(), inplace=True)
df['CO'].fillna(df['CO'].median(), inplace=True)
df['O3'].fillna(df['O3'].median(), inplace=True)
결측 처리 확인
print(round(df.isnull().sum()/len(df.index), 2)*100)
PM2.5 0.0
PM10 0.0
SO2 0.0
NO2 0.0
CO 0.0
O3 0.0
TEMP 0.0
PRES 0.0
DEWP 0.0
RAIN 0.0
wd 0.0
WSPM 0.0
station 0.0
dtype: float64
3. 데이터셋 구성
우선 이번 예제에서는 station을 1개만 선정해 수행하도록 하겠습니다.
station 목록에서 ‘Dongsi’에 대한 내용만 가져온 후 station 칼럼을 삭제하겠습니다.
ex_df = df[df['station']=='Dongsi']
ex_df.drop(['station'],axis=1,inplace=True)
ex_df
| PM2.5 | PM10 | SO2 | NO2 | CO | O3 | TEMP | PRES | DEWP | RAIN | wd | WSPM | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| date | ||||||||||||
| 2013-03-01 00:00:00 | 9.0 | 9.0 | 3.0 | 17.000000 | 300.0 | 89.0 | -0.5 | 1024.5 | -21.4 | 0.0 | NNW | 5.7 |
| 2013-03-01 01:00:00 | 4.0 | 4.0 | 3.0 | 16.000000 | 300.0 | 88.0 | -0.7 | 1025.1 | -22.1 | 0.0 | NW | 3.9 |
| 2013-03-01 02:00:00 | 7.0 | 7.0 | 7.0 | 17.000000 | 300.0 | 60.0 | -1.2 | 1025.3 | -24.6 | 0.0 | NNW | 5.3 |
| 2013-03-01 03:00:00 | 3.0 | 3.0 | 5.0 | 18.000000 | 900.0 | 45.0 | -1.4 | 1026.2 | -25.5 | 0.0 | N | 4.9 |
| 2013-03-01 04:00:00 | 3.0 | 3.0 | 7.0 | 50.638586 | 200.0 | 84.0 | -1.9 | 1027.1 | -24.5 | 0.0 | NNW | 3.2 |
| … | … | … | … | … | … | … | … | … | … | … | … | … |
| 2017-02-28 19:00:00 | 16.0 | 51.0 | 3.0 | 29.000000 | 400.0 | 73.0 | 12.5 | 1013.5 | -16.2 | 0.0 | NW | 2.4 |
| 2017-02-28 20:00:00 | 18.0 | 45.0 | 3.0 | 43.000000 | 500.0 | 54.0 | 11.6 | 1013.6 | -15.1 | 0.0 | WNW | 0.9 |
| 2017-02-28 21:00:00 | 23.0 | 58.0 | 5.0 | 61.000000 | 700.0 | 28.0 | 10.8 | 1014.2 | -13.3 | 0.0 | NW | 1.1 |
| 2017-02-28 22:00:00 | 23.0 | 53.0 | 9.0 | 75.000000 | 900.0 | 15.0 | 10.5 | 1014.4 | -12.9 | 0.0 | NNW | 1.2 |
| 2017-02-28 23:00:00 | 30.0 | 71.0 | 11.0 | 87.000000 | 1200.0 | 4.0 | 8.6 | 1014.1 | -15.9 | 0.0 | NNE | 1.3 |
35064 rows × 12 columns
칼럼 중 wd에 대한 내용은 문자로 되어있어 sklearn에서 제공하는 LabelEncoder을 사용해 인코딩해줍니다.
encoder = LabelEncoder()
ex_df['wd'] = encoder.fit_transform(ex_df['wd'])
| PM2.5 | PM10 | SO2 | NO2 | CO | O3 | TEMP | PRES | DEWP | RAIN | wd | WSPM | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| date | ||||||||||||
| 2013-03-01 00:00:00 | 9.0 | 9.0 | 3.0 | 17.000000 | 300.0 | 89.0 | -0.5 | 1024.5 | -21.4 | 0.0 | 6 | 5.7 |
| 2013-03-01 01:00:00 | 4.0 | 4.0 | 3.0 | 16.000000 | 300.0 | 88.0 | -0.7 | 1025.1 | -22.1 | 0.0 | 7 | 3.9 |
| 2013-03-01 02:00:00 | 7.0 | 7.0 | 7.0 | 17.000000 | 300.0 | 60.0 | -1.2 | 1025.3 | -24.6 | 0.0 | 6 | 5.3 |
| 2013-03-01 03:00:00 | 3.0 | 3.0 | 5.0 | 18.000000 | 900.0 | 45.0 | -1.4 | 1026.2 | -25.5 | 0.0 | 3 | 4.9 |
| 2013-03-01 04:00:00 | 3.0 | 3.0 | 7.0 | 50.638586 | 200.0 | 84.0 | -1.9 | 1027.1 | -24.5 | 0.0 | 6 | 3.2 |
| … | … | … | … | … | … | … | … | … | … | … | … | … |
| 2017-02-28 19:00:00 | 16.0 | 51.0 | 3.0 | 29.000000 | 400.0 | 73.0 | 12.5 | 1013.5 | -16.2 | 0.0 | 7 | 2.4 |
| 2017-02-28 20:00:00 | 18.0 | 45.0 | 3.0 | 43.000000 | 500.0 | 54.0 | 11.6 | 1013.6 | -15.1 | 0.0 | 14 | 0.9 |
| 2017-02-28 21:00:00 | 23.0 | 58.0 | 5.0 | 61.000000 | 700.0 | 28.0 | 10.8 | 1014.2 | -13.3 | 0.0 | 7 | 1.1 |
| 2017-02-28 22:00:00 | 23.0 | 53.0 | 9.0 | 75.000000 | 900.0 | 15.0 | 10.5 | 1014.4 | -12.9 | 0.0 | 6 | 1.2 |
| 2017-02-28 23:00:00 | 30.0 | 71.0 | 11.0 | 87.000000 | 1200.0 | 4.0 | 8.6 | 1014.1 | -15.9 | 0.0 | 5 | 1.3 |
35064 rows × 12 columns
그 후, 학습 수행을 위한 각 칼럼 값들에 대해 0~1사이의 값으로 스케일링을 진행해줍니다.
최종 과정에서의 편의를 위해 라벨 스케일러와 데이터 스케일러를 나누어 사용했습니다.
필요에 따라 하나로 사용할 경우, 이후 원래 범위대로 돌리는 과정에서 데이터 칼럼의 순서를 잘 맞추어 진행하면 됩니다.
# label
label_scaler = MinMaxScaler(feature_range=(0,1))
scaled_label = label_scaler.fit_transform(ex_df['PM2.5'].values.reshape(-1,1))
# data
data_scaler = MinMaxScaler(feature_range=(0,1))
data_df = ex_df.copy()
data_df.drop('PM2.5', axis=1, inplace=True)
scaled_data = data_scaler.fit_transform(data_df)
# concated data(label+data)
scaled = np.concatenate((scaled_label, scaled_data), axis=1)
이렇게 변환된 데이터를 가지고, 시계열 데이터 입력에 맞추어 변형을 진행합니다.
해당 과정에서는 예측 수행 시 살펴볼 시점과, 예측할 시점에 따라 데이터를 재구성하는 것입니다.
즉, 만일 내가 4시간 전까지의 데이터를 모두 사용해 현재 시점을 예측하고자 하는경우, t-4, t-3, t-2, t-1, t 이런식으로 각 칼럼의 해당 시점들의 정보들을 가지고 재구성하는 것 입니다.
해당 내용에 대한 자세한 설명은 다음 링크를 참고하세요.
https://diane-space.tistory.com/316
def series2Supervised(data, n_in=1, n_out=1, dropnan=True):
n_vars = 1 if type(data) is list else data.shape[1]
df = pd.DataFrame(data)
cols, names = list(), list()
# 입력값 순서 (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)]
# 예측값 순서(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)]
# 합치기
aggregate = pd.concat(cols, axis=1)
aggregate.columns=names
# nan 값 제거
if dropnan:
aggregate.dropna(inplace=True)
return aggregate
우리는 바로 직전 시점(t-1)을 사용해 t를 예측하도록 할 예정입니다.
reframed = series2Supervised(scaled,1, 1)
reframed
| var1(t-1) | var2(t-1) | var3(t-1) | var4(t-1) | var5(t-1) | var6(t-1) | var7(t-1) | var8(t-1) | var9(t-1) | var10(t-1) | … | var3(t) | var4(t) | var5(t) | var6(t) | var7(t) | var8(t) | var9(t) | var10(t) | var11(t) | var12(t) | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0.008174 | 0.007345 | 0.009057 | 0.058594 | 0.020202 | 0.082549 | 0.281520 | 0.681239 | 0.216849 | 0.0 | … | 0.009057 | 0.054688 | 0.020202 | 0.081615 | 0.278066 | 0.692168 | 0.205928 | 0.0 | 0.4375 | 0.371429 |
| 2 | 0.001362 | 0.002099 | 0.009057 | 0.054688 | 0.020202 | 0.081615 | 0.278066 | 0.692168 | 0.205928 | 0.0 | … | 0.022403 | 0.058594 | 0.020202 | 0.055456 | 0.269430 | 0.695811 | 0.166927 | 0.0 | 0.3750 | 0.504762 |
| 3 | 0.005450 | 0.005247 | 0.022403 | 0.058594 | 0.020202 | 0.055456 | 0.269430 | 0.695811 | 0.166927 | 0.0 | … | 0.015730 | 0.062500 | 0.080808 | 0.041442 | 0.265976 | 0.712204 | 0.152886 | 0.0 | 0.1875 | 0.466667 |
| 4 | 0.000000 | 0.001049 | 0.015730 | 0.062500 | 0.080808 | 0.041442 | 0.265976 | 0.712204 | 0.152886 | 0.0 | … | 0.022403 | 0.189994 | 0.010101 | 0.077878 | 0.257340 | 0.728597 | 0.168487 | 0.0 | 0.3750 | 0.304762 |
| 5 | 0.000000 | 0.001049 | 0.022403 | 0.189994 | 0.010101 | 0.077878 | 0.257340 | 0.728597 | 0.168487 | 0.0 | … | 0.029076 | 0.089844 | 0.020202 | 0.072272 | 0.248705 | 0.735883 | 0.218409 | 0.0 | 0.4375 | 0.228571 |
| … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … |
| 35059 | 0.013624 | 0.047219 | 0.005720 | 0.062500 | 0.020202 | 0.079747 | 0.521589 | 0.471767 | 0.308892 | 0.0 | … | 0.009057 | 0.105469 | 0.030303 | 0.067601 | 0.506045 | 0.480874 | 0.297972 | 0.0 | 0.4375 | 0.228571 |
| 35060 | 0.017711 | 0.051417 | 0.009057 | 0.105469 | 0.030303 | 0.067601 | 0.506045 | 0.480874 | 0.297972 | 0.0 | … | 0.009057 | 0.160156 | 0.040404 | 0.049850 | 0.490501 | 0.482696 | 0.315133 | 0.0 | 0.8750 | 0.085714 |
| 35061 | 0.020436 | 0.045121 | 0.009057 | 0.160156 | 0.040404 | 0.049850 | 0.490501 | 0.482696 | 0.315133 | 0.0 | … | 0.015730 | 0.230469 | 0.060606 | 0.025559 | 0.476684 | 0.493625 | 0.343214 | 0.0 | 0.4375 | 0.104762 |
| 35062 | 0.027248 | 0.058762 | 0.015730 | 0.230469 | 0.060606 | 0.025559 | 0.476684 | 0.493625 | 0.343214 | 0.0 | … | 0.029076 | 0.285156 | 0.080808 | 0.013414 | 0.471503 | 0.497268 | 0.349454 | 0.0 | 0.3750 | 0.114286 |
| 35063 | 0.027248 | 0.053515 | 0.029076 | 0.285156 | 0.080808 | 0.013414 | 0.471503 | 0.497268 | 0.349454 | 0.0 | … | 0.035749 | 0.332031 | 0.111111 | 0.003137 | 0.438687 | 0.491803 | 0.302652 | 0.0 | 0.3125 | 0.123810 |
35035 rows × 24 columns
해당 칼럼에서 우리는 PM2.5를 타겟으로 사용할 예정입니다. 이외의 t시점에서의 다른 칼럼들은 필요없으므로 제거해 줍니다.
reframed.drop(reframed.columns[-11:], axis=1, inplace=True)
학습 데이터셋 / validation 데이터셋 / 테스트 데이터셋 총 3개로 구성해 진행합니다.
우리는 학습을 위해 3년치 데이터를 사용할 예정입니다.
values = reframed.values
n_train_hours = (24*30*11+24*28)*3
train = values[:n_train_hours, :]
test = values[n_train_hours:, :]
train_x, train_y = train[:, :-1], train[:, -1]
test_x, test_y = test[:, :-1], test[:, -1]
# reshape
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(f'train:{train_x.shape}/{train_y.shape}\ntest:{test_x.shape}/{test_y.shape}')
# evaluation dataset 생성
x_train, x_valid, y_train, y_valid = train_test_split(train_x, train_y, test_size=0.2, shuffle=False)
train:(25776, 1, 12)/(25776,)
test:(9259, 1, 12)/(9259,)
4. 학습
모델 구성 LSTM 레이어 1개를 가지고 수행할 예정입니다.
모델 구성의 경우 하이퍼파라미터를 다양하게 테스트해 최적의 성능으로 개선가능합니다.
이번 모델의 경우, unit=128, loss=mae, optimizer=adam, activation=tanh로 구성했습니다.
model = tf.keras.models.Sequential()
model.add(tf.keras.layers.Input((train_x.shape[1], train_x.shape[2])))
model.add(tf.keras.layers.LSTM(128, activation='tanh'))
model.add(tf.keras.layers.Dense(1))
model.compile(loss='mae', optimizer='adam')
model.summary()
Model: "sequential"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
lstm (LSTM) (None, 128) 72192
_________________________________________________________________
dense (Dense) (None, 1) 129
=================================================================
Total params: 72,321
Trainable params: 72,321
Non-trainable params: 0
_________________________________________________________________
과적합 방지를 위해 early_stopping을 적용해 학습을 수행한 결과 epoch를 40회 채우지 못한 채 종료되었습니다.
early_stop = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=4)
history = model.fit(x_train, y_train, epochs=40, batch_size=32,
validation_data=(x_valid, y_valid),
verbose=2,
shuffle=False,
callbacks=[early_stop])
5. 테스트
이제 확인을 위해 예측을 수행하겠습니다.
pred = model.predict(test_x)
원본과의 확인을 위해 정규화 전 데이터로 데이터를 역변환 합니다.
해당 내용은 이후 성능 평가에서 원본과의 비교를 위해 진행되며, 값의 흐름만을 비교하고자 할때는 생략해도 됩니다.
test_x = test_x.reshape((test_x.shape[0], test_x.shape[2])) # 원래 모양으로 되돌리기
inv_pred = label_scaler.inverse_transform(pred)
test_y = test_y.reshape((len(test_y),1))
inv_y = label_scaler.inverse_transform(test_y)
plt.figure(figsize=(30,10))
plt.plot(inv_y, label='label')
plt.plot(inv_pred, label='predict')
plt.legend()
plt.show()
모델 성능 평가 결과입니다.
RMSE: 21.786
6. 추가 성능 비교
추가적으로 LSTM 레이어를 여러 겹 중첩한 모델과의 비교를 진행했습니다.
여러 겹 중첩하는 경우, 다른 파라미터는 고정시킨 채 레이어 층 수만을 추가해 비교 진행했습니다.
(모델 생략)
| Layer 수 | Train loss / Validation loss | RMSE | 비고 |
|---|---|---|---|
| 2 | 0.0178 / 0.0167 | 22.970 | epoch 16 |
| 3 | 0.0206 / 0.0238 | 29.024 | epoch 8 |
| 4 | 0.0212 / 0.0254 | 29.735 | epoch 7 |
테스트 결과, 중첩을 진행한 경우, 레이어 2개일 때가 가장 성능이 좋게 나왔으며, 중첩하지 않은 테스트 결과가 가장 좋아 무조건적으로 중첩을 수행한다고 성능이 나아지지 않음을 보였습니다.
본 포스팅에서는 LSTM을 사용해 베이징의 공기질 예측을 수행한 예제를 수행했습니다.
데이터셋을 불러오고, 결측치를 채운 후, 시계열 데이터셋으로 재구성하는 과정에 대해 진행했습니다.
가장 기본적인 LSTM 레이어 1개로 구성된 모델과 추가적으로 Layer을 중첩한 모델까지 구성해 테스트를 진행했습니다.
타겟 데이터인 PM2.5는 데이터의 범위가 3~681임을 감안했을 때, rmse값이 21.786으로 엄청 예측을 못한다고는 볼 수 없었습니다.
다양하게 모델을 구성해 하이퍼파라미터를 찾아간다면, 더 나은 성능을 보일 수 있을 것입니다.
이상으로 본 포스팅을 마무리하겠습니다.