Help us understand the problem. What is going on with this article?

ディープラーニングで7月の気温を分析(3) -CNNで回帰-

More than 1 year has passed since last update.

前回までは、過去の最高気温から、次の日の温度クラスを推定する分類問題を
扱ってきました。今回は趣向を変えて、CNNを使って最高気温を直接推定する
回帰問題を扱いたいと思います。

figure_11.png

問題設定

問題設定を再定義します。以下の図のような問題を考えます。
温度.png

過去5日分の気温推移から次の日の最高気温を予測するものです。

今回の出力は回帰ですので、「7/6は30℃」みたいな形で出てきます。
ただし、気温は天気の影響を大きく受けるため、yを正確に予測するには
その日の天気を正確に予測しないといけません。

しかし、過去5日分の気温推移から、次の日の天気を予測することはほぼ不可能です。
そこで、出力yは「予測値は30℃だけど、天気によって±5℃くらい変動するよ」と
いう設定に置き換えます。

そして、30±5℃から外れたものを異常、つまり特異な温度推移をしている
と見なして異常検知させます。下の図を例にします。

figure_8.png

青い線は実測値で、赤い実線は予測値です。点線が変動範囲を示しています。
図中の赤い矢印は変動範囲内なので正常ですが、青い矢印は範囲外なので
異常と見なします。今回こそ、2018年の高温が異常と出れば良いのですが。。。

今回の問題は、次の日の最高気温をT、予測値をt、変動幅をeとすると
以下の式になります。

T=t±e

そして、Tから外れたものを異常と見なします。

このような問題を解く場合、ベイズ線形回帰問題↓として扱えば、きれいにかつ
早くできるのでしょうが、今回はCNNを使います。
https://qiita.com/naoya_t/items/80ea108cebc694f5cd63

CNNモデルの構築

今回の問題を解くのにあたり、以下の図の方式でCNNを学習させます。

CNNベイズ線形回帰.png

・まずは、過去5日分の気温推移でCNN1を学習させ、気温を直接推定させます。(図中の黄色)
・そして、CNN1の予測値と実測値の差をとります。(図中の緑)
・最後に、過去5日分の気温推移でCNN2を学習させ、偏差eを推定させます。(図中のオレンジ)

ただし、正常と見なす基準は図中のeを使って以下のように設定します。

T=t±1.5e

ちょっと基準を広くしました。

kerasで実装

まずは、データを取得して、時系列に並び替えます。

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np

#30年分の 温度データロード,整列
X = np.loadtxt("temperature.csv",delimiter=",")
data = X[:20]
for i in range(1,30):
    x = X[i*20:i*20+20]
    data = np.vstack((data,x))

#テストデータと学習データの分割
x_train0 = data[:20]
x_test0 = data[20:]

from sklearn.cross_validation import train_test_split
#新バージョンはfrom sklearn.cross_model_selection

#時系列データの作成
def make_data(x_train_rnn,Lookback=5):
    data, target = [], []

    for i in range(x_train_rnn.shape[0]):
        for j in range(x_train_rnn.shape[1]-Lookback):
            data.append(x_train_rnn[i,j:j+Lookback])
            target.append(x_train_rnn[i,j+Lookback])

    re_data = np.array(data).reshape(len(data), Lookback, 1)
    re_target = target

    return re_data, re_target

x_train, y_train = make_data(x_train0)
X_test, Y_test = make_data(x_test0)

#trainとvalidationを分割
X_train, X_val, Y_train, Y_val = train_test_split(x_train, y_train, test_size=0.2)

#入力が0~1になるように修正
x_max = np.max(X_train)
x_min = np.min(X_train)

X_train = (X_train-x_min)/(x_max-x_min)
Y_train = (Y_train-x_min)/(x_max-x_min)
X_val = (X_val-x_min)/(x_max-x_min)
Y_val = (Y_val-x_min)/(x_max-x_min)
X_test = (X_test-x_min)/(x_max-x_min)
Y_test = (Y_test-x_min)/(x_max-x_min)

#データを4次元化
lookback = 5
X_train = X_train.reshape((len(X_train),lookback,1,1))
X_val = X_val.reshape((len(X_val),lookback,1,1))
X_test = X_test.reshape((len(X_test),lookback,1,1))

次に、CNN1を構築します。

from keras.models import Model
from keras.layers import Dense, Dropout, Activation, Flatten, Add, Input
from keras.layers.advanced_activations import LeakyReLU
from keras.layers.convolutional import Conv2D, MaxPooling2D
from keras.layers.normalization import BatchNormalization

#resnet
def resblock(x, filters, kernel_size):
    x_ = Conv2D(filters, kernel_size, padding='same')(x)
    x_ = BatchNormalization()(x_)
    x_ = Activation(LeakyReLU())(x_)
    x_ = Conv2D(filters, kernel_size, padding='same')(x_)
    x = Add()([x_, x])
    x = BatchNormalization()(x)
    x = Activation(LeakyReLU())(x)
    return x

#CNN1(気温の予測値)の学習
input_ = Input(shape=(lookback, 1,1))#横の数、縦の数、RGB

c = Conv2D(64, (2, 1),padding='same')(input_)
x = BatchNormalization()(c)
c = Activation(LeakyReLU())(c)
c = Dropout(0.2)(c)
c = resblock(c,filters=64, kernel_size=(2, 1))
c = MaxPooling2D(pool_size=(2, 1))(c)

c = Flatten()(c)
c = Dense(30)(c)
c = Activation(LeakyReLU())(c)
c = Dropout(0.2)(c)
c = Dense(1, activation='sigmoid')(c)

model = Model(input_, c)

model.compile(loss='mse', optimizer='adam', metrics=['mse'])

hist = model.fit(X_train, Y_train, batch_size = 10, epochs=50, verbose=1, shuffle=True,
                 validation_data = (X_val,Y_val))

#学習結果描画
plt.figure()               
plt.plot(hist.history['val_loss'],label="val_loss")
plt.plot(hist.history['loss'],label="train_loss")
plt.legend()
plt.show()

損失関数の推移↓
figure_1.jpg

次に、CNN2を構築します。

#偏差の計算
predict = model.predict(X_train,batch_size=1)
predict = predict.reshape(len(predict))
Y_train_e = np.absolute(Y_train - predict)

predict_val = model.predict(X_val,batch_size=1)
predict_val = predict_val.reshape(len(predict_val))
Y_val_e = np.absolute(Y_val - predict_val)

#入力が0~1になるように修正
e_max = np.max(Y_train_e)
e_min = np.min(Y_train_e)

Y_train_e = (Y_train_e-e_min)/(e_max - e_min)
Y_val_e = (Y_val_e-e_min)/(e_max - e_min)

#CNN2(編差)の学習
input_ = Input(shape=(lookback, 1,1))#横の数、縦の数、RGB

c = Conv2D(64, (2, 1),padding='same')(input_)
c = BatchNormalization()(c)
c = Activation(LeakyReLU())(c)
c = Dropout(0.2)(c)
c = resblock(c,filters=64, kernel_size=(2, 1))
c = MaxPooling2D(pool_size=(2, 1))(c)

c = Flatten()(c)
c = Dense(30)(c)
c = Activation(LeakyReLU())(c)
c = Dropout(0.2)(c)
c = Dense(1, activation='sigmoid')(c)

model_e = Model(input_, c)

model_e.compile(loss='mse', optimizer='adam', metrics=['mse'])

hist = model_e.fit(X_train, Y_train_e, batch_size = 10, epochs=50, verbose=1, shuffle=True,
                 validation_data = (X_val,Y_val_e))

#学習結果描画
plt.figure()               
plt.plot(hist.history['val_loss'],label="val_loss")
plt.plot(hist.history['loss'],label="train_loss")
plt.legend()
plt.show()

損失関数の推移↓
figure_2.jpg

結果と考察

各年の実測値と予測値を描画します。
特徴的なグラフを見ていきます。

#予測範囲描画
t_predict = model.predict(X_test,batch_size=1)
t_predict = (x_max-x_min)*t_predict + x_min
t_predict = t_predict.reshape((x_test0.shape[0],x_test0.shape[1]-5,))

t_predict_e = model_e.predict(X_test,batch_size=1)
t_predict_e = 1.5*((x_max-x_min)*((e_max-e_min)*t_predict_e+e_min))
t_predict_e = t_predict_e.reshape((x_test0.shape[0],x_test0.shape[1]-5,))

for i in range(len(x_test0)):
    plt.figure()
    year = 2009+i
    plt.plot(range(1,21),x_test0[i], label="true")
    plt.plot(range(6,21),t_predict[i],c="r",label="predict")
    plt.plot(range(6,21),t_predict[i] + t_predict_e[i], c="r",linestyle="dashed")
    plt.plot(range(6,21),t_predict[i] - t_predict_e[i], c="r", linestyle="dashed")
    plt.title(year)
    plt.ylim(20,40)
    plt.xlim(1,20)
    plt.xlabel("date")
    plt.ylabel("temperature")
    plt.legend(loc="lower right")
    plt.show()

figure_11.png

図の縦軸は温度(℃)、横軸は日付(左端は7/1、右端は7/20)です。
2017年は、予測値がほとんど一致しており、例年通りの推移だったと
予想されます。

figure_4.png

2010年は、7/13の低温は異常と認識しています。さらに、7/17の右肩
上がりの気温上昇も異常と認識しています。

figure_10.png

2016年は、7/9の気温急降下は異常でしたが、7/11は予測値と実測値が
ピタリと一致し、かつ変動幅が小さなっています。CNNは余程自信が
あったのでしょう。

figure_9.png

2015年は、7/6~9の低温期が異常と認識できるのかが焦点でした。
そして、その期間はほとんど異常と認識しています。

figure_12.png

2018年は、7/15以降の高温期が異常と認識できるのかが焦点でしたが
異常と認識できたのは、7/15と7/20の二日だけでした。
惜しい!といったところでしょうか。

一応、全テストデータの異常点を赤いプロットで示してみます。

#異常推移を赤いプロットで描画
plot_anomaly_x = []
plot_anomaly_y = []

for i in range(x_test0.shape[0]):
    for j in range(x_test0.shape[1]):
        if j > 4:#6日以降の予測値で判定
            if t_predict[i,j-5] + t_predict_e[i,j-5] < x_test0[i,j] or t_predict[i,j-5] - t_predict_e[i,j-5] > x_test0[i,j]:
                plot_anomaly_x.append(j)
                plot_anomaly_y.append(x_test0[i,j])

plt.figure()
for i in range(len(x_test0)):
    year = 2009+i
    plt.plot(range(1,21),x_test0[i], color=cm.jet(float(i) / len(x_test0)),label=year)
#plt.legend(loc="lower right")
plt.scatter(np.array(plot_anomaly_x)+1,np.array(plot_anomaly_y),c="r")
plt.title("Test data")
plt.ylim(20,40)
plt.xlim(5,20)
plt.xlabel("date")
plt.ylabel("temperature")
plt.show()

figure_13.png

前回より、異常と認識している点が増えてしまいましたが、納得いく
異常点が増えたかなぁと思います。

まとめ

これまで3回にわたり、ディープラーニングを使って7月の気温を
分析してきました。以下に、感想をまとめました。

・時系列データに対しては、LSTMよりCNNの方が若干優れているような
 気がします。

・異常検知させるに当たり、分類と回帰のどちらを用いるべきかという
 問いに関しては、回帰が良い気がします。

 回帰をお勧めする1つ目の理由は、設定の簡単さです。分類だとクラス分け
 作業が発生するため、クラス分けの基準を決めないといけません。その作業が
 面倒です。さらに、基準を都合の良いように動かすことも可能なので、
 解析結果を本当に信じて良いか分からなくなってしまいます。
 その点、回帰は予測値がドン!と出てくるので、簡単かつ都合よく
 動かすことはできません。

 回帰をお勧めする2つ目の理由は、直感的に分かりやすいという点です。
 分類は確率しか出てこないので、その数値を信じるしかありませんが、
 回帰は予想値がそのまま出てくるので扱いやすいです。さらに、
 回帰+偏差を用いるとCNNの予測した幅が可視化できて分かりやすいです。

 時系列データに対し、「確率が・・・」と説明しても客先や上司には、なかなか
 伝わりにくいです。しかし、「回帰で予測した幅がこれで、実測値がその中に
 入っているでしょ」というと圧倒的に分かりやすいです。このアドバイスは
 cranpunさんから頂きました。ありがとうございます。

 以上、長々と書いてしまいましたが、論文調査をしっかり行っていない
 ので、話し半分に聞いて下さい。

shinmura0
自己紹介はツイッターをご覧ください。 https://twitter.com/shinmura0
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away