RNN/LSTMで気温予測【3. 多変数編】

はじめに

前回は気温のみのデータから予測を行うモデルを構築した。

今回は他のデータも用いることで精度の向上を図る。

コードの変更

Datasetの次元が変わることになる。n_featuresを特徴量の数としてところどころコードを変更する。

Dataset/DataLoaderの作成

from sklearn.preprocessing import MinMaxScaler

import torch
from torch.utils.data import DataLoader

# 使用する特徴量
features = ['temperature', 'atmosphere1', 'dew point', 'humidity']
n_features = len(features)

y = df[features].values
x = np.arange(len(y))

train_rate = 0.95
train_size = int(len(y)*train_rate)
test_size = len(y) - train_size

# 最近の傾向を学習するため訓練データを新しい95%分、テストデータを古い5%分に変更している。
y_train = y[-train_size:]
y_test = y[:test_size]

print("train size: {}, test size: {} ".format(len(y_train), len(y_test)))

ms = MinMaxScaler()
y = y.reshape(-1, n_features)
y_train = y_train.reshape(-1, n_features)
y_test = y_test.reshape(-1, n_features)
y_ms = ms.fit_transform(y)
y_train_ms = ms.transform(y_train)
y_test_ms = ms.transform(y_test)

time_step = 120
BATCH_SIZE = 1024
BATCH_SIZE_test = 512

n_sample = len(y) - time_step
n_sample_train = train_size - time_step
n_sample_test = test_size - time_step

input_data = np.zeros((n_sample, time_step, n_features))
input_data_train = np.zeros((n_sample_train, time_step, n_features))
correct_data_train = np.zeros((n_sample_train, n_features))
input_data_test = np.zeros((n_sample_test, time_step, n_features))
correct_data_test = np.zeros((n_sample_test, n_features))

for i in range(n_sample):
    input_data[i] = y_ms[i:i+time_step]

for i in range(n_sample_train):
    input_data_train[i] = y_train_ms[i:i+time_step]
    correct_data_train[i] = y_train_ms[i+time_step]

for i in range(n_sample_test):
    input_data_test[i] = y_test_ms[i:i+time_step]
    correct_data_test[i] = y_test_ms[i+time_step]
    
input_data_tensor = torch.tensor(input_data, dtype=torch.float)
input_data_train_tensor = torch.tensor(input_data_train, dtype=torch.float)
correct_data_train_tensor = torch.tensor(correct_data_train, dtype=torch.float)
input_data_test_tensor = torch.tensor(input_data_test, dtype=torch.float)
correct_data_test_tensor = torch.tensor(correct_data_test, dtype=torch.float)

dataset_train = torch.utils.data.TensorDataset(input_data_train_tensor, correct_data_train_tensor)
dataset_test = torch.utils.data.TensorDataset(input_data_test_tensor, correct_data_test_tensor)
train_loader = DataLoader(dataset_train, batch_size=BATCH_SIZE, shuffle=True)
test_loader = DataLoader(dataset_test, batch_size=BATCH_SIZE_test, shuffle=True)

予測

# 何時間後まで予測するか
pred_range = 24

for i in range(3):
    # 予測開始時刻、訓練データ内で指定
    hour_id = np.random.randint(time_step, 240000)
    predicted = input_data_train_tensor[hour_id].numpy()
    model.eval()
    with torch.no_grad():
        for k in range(pred_range):
            x = torch.tensor(predicted[-time_step:], dtype=torch.float)
            x = x.reshape(1, time_step, n_features)
            pred_y = model(x)
            predicted = np.append(predicted, pred_y.cpu().numpy(), axis=0)
    predicted = predicted.reshape(-1, n_features)
    predicted = ms.inverse_transform(predicted)

    fig = plt.figure(figsize=(20, 10))
    plt.subplots_adjust(wspace=0.2, hspace=0.2)
    for j in range(n_features):
        ax = fig.add_subplot(2, 2, j+1)
        ax.plot(range(test_size+hour_id, test_size+hour_id+time_step+pred_range), y[(test_size+hour_id):(test_size+hour_id+time_step+pred_range), j], label='Correct')
        ax.plot(range(test_size+hour_id, test_size+hour_id+time_step+pred_range), predicted[:, j], label='pred_y')
        ax.set_title(features[j])
        ax.legend()
    filename = 'model' + str(model_num) + '/pred_sample_' + str(i) + '.jpg' 
    plt.savefig(filename, dpi=100)
    plt.show()

結果の考察用

div = 5
pred_hour = 24
n_data = int((n_sample_test - pred_hour) / div) + 1
model.eval()
loss_data = np.zeros((n_features, pred_hour, 1))

for i in range(n_data):
    print('\r i: {}'.format(i*div), end='')
    predicted = input_data_test_tensor[div*i].numpy()
    with torch.no_grad():
        for j in range(pred_hour):
            x = torch.tensor(predicted[-time_step:], dtype=torch.float)
            x = x.reshape(1, time_step, n_features)
            pred_y = model(x)
            predicted = np.append(predicted, pred_y.cpu(), axis=0)
    predicted = predicted.reshape(-1, n_features)
    predicted = ms.inverse_transform(predicted)
    loss_ = abs(predicted - y_test[div*i:div*i+time_step+pred_hour])
    loss_data = np.append(loss_data, loss_[-pred_hour:].T.reshape(n_features, pred_hour, 1), axis=2)

loss_data = np.transpose(loss_data[:, :, 1:], (0, 2, 1))
zero_series = np.zeros((n_features, n_data, 1))
loss_data = np.append(zero_series, loss_data, axis=2)
loss_df = pd.DataFrame(loss_data[0])

パラメータ

パラメータのうち以下のコードを変更する

n_inputs = n_features
n_outputs = n_features

結果

model1

特徴量は気温、気圧、湿度とした。

features = ['temperature', 'atmosphere1', 'humidity']
# Hyperparameters -------------------------------------
time_step = 120
BATCH_SIZE = 1024
BATCH_SIZE_test = 512
n_hidden = 64
n_layers = 1
LEARNING_RATE = 0.00005
EPOCHS = 300 
# -----------------------------------------------------

損失関数の値の推移

予測の様子

気温や湿度に関しては比較的うまく予測できているが、気圧の予測は大きく外れているように見える天気に関わりがある指標なのでこれがもう少し上手くいけば気温予測もよりよくできるかもしれない。天気は西から変化することが多いので東京以外の都市のデータを用いれば精度が上がるかもしれない。

時間ごとの誤差に関する統計

24時間後でも平均2度以内を達成した。中央値や75%の値も1変数のときより低下している。

model2

今度は露点温度のデータも追加した。また、schedulerを用いて途中で学習率を減衰させることで精度の向上をねらった。

features = ['temperature', 'atmosphere1', 'dew point', 'humidity']
# Hyperparameters -------------------------------------
time_step = 240
BATCH_SIZE = 1024
BATCH_SIZE_test = 512
n_hidden = 512
n_layers = 1
LEARNING_RATE = 0.01
EPOCHS = 500
# -----------------------------------------------------

scheduler = optim.lr_scheduler.MultiStepLR(optimizer, milestones=[50, 100, 150, 200, 250, 300, 350, 400, 450], gamma=0.2)

エポックのループにscheduler.step()も記述する。

損失関数の値の推移

縮尺の問題で見づらいが、スケジューラのおかげで50エポックや100エポックの部分で急激に訓練データのlossが減少し、訓練データの振動が穏やかになった。一方で300以降でのマイルストーンでは大きな変化は見られず、ともに0.000543程度が収束値となった。

予測の様子

時間ごとの誤差に関する統計

24時間後の予測で平均1.77度、75%は2.46度以内に収まった。model1よりも精度が向上したといえるだろう。

終わりに

次回は他の地点のデータを利用することでより高い精度を出せるようなモデルを作成したい。