IDWR データで再帰型ニューラルネットワーク - Keras

Python 時系列データ DeepLearning

前回 加工した IDWR データ を使って再帰型ニューラルネットワーク（RNN）を Keras + TensorFlow で試してみました。

ソースは http://github.com/fits/try_samples/tree/master/blog/20180121/

はじめに

インストール

tensorflow と keras をインストールしておきます。

インストール例
> pip install tensorflow
> pip install keras

データセット

データセットは 前回 に作成した、2014年 1週目 ～ 2017年 49週目 （2015年は 53週目まで） の teiten.csv を 1つの csv ファイル （idwr.csv） へ連結したものを使います。

再帰型ニューラルネットワーク（RNN）

指定の感染症に対する週別の全国合計値を再帰型ニューラルネットワーク（RNN）で学習・予測する処理を実装してみます。

(a) 4週分のデータを使用

とりあえず、任意の週の過去 4週の報告数を入力データ、その週の報告数をラベルデータとして学習を実施してみます。

ただ、入力データに季節性（周期性）の指標となるデータ（ここでは週）を含んでいないので、季節性のあるデータに対する結果は期待できないように思います。

データ加工

入力データとラベルデータは、以下の要領で 1週分ずつスライドしたものを用意します。

  • 2014年 1 ～ 4週目を入力データとし 5週目をラベルデータ
  • 2014年 2 ～ 5週目を入力データとし 6週目をラベルデータ
入力データ・ラベルデータの作成例（インフルエンザ）
import pandas as pd
import numpy as np

t = 4

df = pd.read_csv('idwr.csv', encoding = 'UTF-8')

ds = df.groupby(['year', 'week'])['インフルエンザ'].sum()

def window(d, wsize):
    return [d[i:(i + wsize)].values.flatten() for i in range(len(d) - wsize + 1)]

# t + 1 週分のスライドデータを作成
dw = window(ds.astype('float'), t + 1)

# 入力データ（t週分）
data = np.array([i[0:-1] for i in dw]).reshape(len(dw), t, 1)
# ラベルデータ（残り 1週分）
labels = np.array([i[-1] for i in dw]).reshape(len(dw), 1)

上記では 1週分ずつスライドさせた 5週分の配列を作成し、4週分を入力データ、残りの 1週分をラベルデータとしています。

その結果、data・labels 変数の内容は以下のようになります。

data 変数の内容
array([[[  9.89100000e+03],
        [  2.71000000e+04],
        [  5.82330000e+04],
        [  1.22618000e+05]],

       [[  2.71000000e+04],
        [  5.82330000e+04],
        [  1.22618000e+05],
        [  1.70403000e+05]],

       [[  5.82330000e+04],
        [  1.22618000e+05],
        [  1.70403000e+05],
        [  1.51829000e+05]],

       ・・・

       [[  2.40700000e+03],
        [  2.58800000e+03],
        [  3.79900000e+03],
        [  7.28000000e+03]],

       [[  2.58800000e+03],
        [  3.79900000e+03],
        [  7.28000000e+03],
        [  1.27850000e+04]]])
labels 変数の内容
array([[  1.70403000e+05],
       [  1.51829000e+05],
       [  1.39162000e+05],
       ・・・
       [  1.27850000e+04],
       [  2.01270000e+04]])

学習

上記のデータを使って再帰型ニューラルネットワークの学習処理を実施してみます。

今回は GRU（ゲート付き再帰的ユニット）を使いました。

学習処理例（Jupyter Notebook 使用）
%matplotlib inline
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.layers.recurrent import GRU
from keras.optimizers import Nadam

epoch = 3000
batch_size = 52
n_num = 80

model = Sequential()

# 再帰型ニューラルネットワークの GRU 使用
model.add(GRU(n_num, activation = 'relu', input_shape = (t, 1)))

model.add(Dense(1))
model.add(Activation('linear'))

opt = Nadam()

# モデルの構築
model.compile(loss = 'mean_squared_error', optimizer = opt)

# 学習
hist = model.fit(data, labels, epochs = epoch, batch_size = batch_size)

# 誤差の遷移状況をグラフ化
plt.plot(hist.history['loss'])

処理結果は以下の通り。 誤差の遷移を見る限り、学習は進んでいるように見えます。

学習処理結果例
Epoch 1/3000
202/202 [==============================] - 1s 7ms/step - loss: 3604321859.1683
Epoch 2/3000
202/202 [==============================] - 0s 193us/step - loss: 3120768191.3663
・・・
Epoch 2999/3000
202/202 [==============================] - 0s 183us/step - loss: 41551096.5149
Epoch 3000/3000
202/202 [==============================] - 0s 193us/step - loss: 41367050.6931

f:id:fits:20180121173527p:plain

予測

3つのグラフを描画して学習したモデルの予測能力を比べてみます。

名称 概要
actual 実データ
predict1 実データを入力データにした予測結果（学習時の入力データをそのまま使用）
predict2 実データの先頭のみを使った予測結果

predict2 は以下のように最初だけ実データを使って、予測結果を次の入力データの一部にして予測を繋げていった結果です。

  • 2014年 1 ～ 4週目を入力データとし 5週目を予測
  • 2014年 2 ～ 4週目と 5週目の予測値を入力データとし 6週目を予測
  • 2014年 3 ～ 4週目と 5 ～ 6週目の予測値を入力データとし 7週目を予測
  • 2014年 4週目と 5 ～ 7週目の予測値を入力データとし 8週目を予測
  • 5 ～ 8週目の予測値を入力データとし 9週目を予測

そのため、実際の予測では predict2 の結果が重要となります。

予測処理例（Jupyter Notebook 使用）
from functools import reduce

# 実データの描画
plt.plot(ds.values, label = 'actual')

# 入力データを使った予測（predict1）
res1 = model.predict(data)
# predict1 の結果を描画
plt.plot(range(t, len(res1) + t), res1, label = 'predict1')

# predict2 用の予測処理
def predict(a, b):
    r = model.predict(a[1])

    return (
        np.append(a[0], r), 
        np.append(a[1][:, 1:], np.array([r]), axis = 1)
    )

# 入力データの先頭
fst_data = data[0].reshape(1, t, 1)
# 入力データの先頭を使った予測（predict2）
res2, _ = reduce(predict, range(len(ds) - t), (np.array([]), fst_data))
# predict2 の結果を描画
plt.plot(range(t, len(res2) + t), res2, label = 'predict2')

plt.legend()

インフルエンザのデータに対して適用した結果が以下です。

予測結果（インフルエンザ）

f:id:fits:20180121173550p:plain

predict1 は実データをほぼ再現できていますが、predict2 は全く駄目でした。

やはり、入力データへ週の情報を与えなかった事が原因だと思われます。

備考

上記処理を単一スクリプト化したものが以下です。

sample1.py
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from functools import reduce
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.layers.recurrent import GRU
from keras.optimizers import Nadam

data_file = sys.argv[1]
item_name = sys.argv[2]
dest_file = sys.argv[3]

t = 4
epoch = 3000
batch_size = 52
n_num = 80

df = pd.read_csv(data_file, encoding = 'UTF-8')

ds = df.groupby(['year', 'week'])[item_name].sum()

def window(d, wsize):
    return [d[i:(i + wsize)].values.flatten() for i in range(len(d) - wsize + 1)]

dw = window(ds.astype('float'), t + 1)

data = np.array([i[0:-1] for i in dw]).reshape(len(dw), t, 1)
labels = np.array([i[-1] for i in dw]).reshape(len(dw), 1)

model = Sequential()

model.add(GRU(n_num, activation = 'relu', input_shape = (t, 1)))

model.add(Dense(1))
model.add(Activation('linear'))

opt = Nadam()

model.compile(loss = 'mean_squared_error', optimizer = opt)

hist = model.fit(data, labels, epochs = epoch, batch_size = batch_size)

fig, axes = plt.subplots(1, 2, figsize = (12, 6))

axes[0].plot(hist.history['loss'])

axes[1].plot(ds.values, label = 'actual')

res1 = model.predict(data)

axes[1].plot(range(t, len(res1) + t), res1, label = 'predict1')

def predict(a, b):
    r = model.predict(a[1])

    return (
        np.append(a[0], r), 
        np.append(a[1][:, 1:], np.array([r]), axis = 1)
    )

fst_data = data[0].reshape(1, t, 1)
res2, _ = reduce(predict, range(len(ds) - t), (np.array([]), fst_data))

axes[1].plot(range(t, len(res2) + t), res2, label = 'predict2')

axes[1].legend()

plt.savefig(dest_file)

(b) 年と週を追加したデータを使用

季節性の強いインフルエンザのようなデータに対して (a) の入力データでは無理があった気がするので、次は入力データへ年と週を追加して試してみます。

入力データとラベルデータ

(a) の入力データへラベルデータの年と週を追加します。

  • 2014 と 5、2014年 1 ～ 4週目を入力データとし 5週目をラベルデータ
  • 2014 と 6、2014年 2 ～ 5週目を入力データとし 6週目をラベルデータ
入力データ（data 変数の内容）例
array([[[  2.01400000e+03],
        [  5.00000000e+00],
        [  9.89100000e+03],
        [  2.71000000e+04],
        [  5.82330000e+04],
        [  1.22618000e+05]],

       [[  2.01400000e+03],
        [  6.00000000e+00],
        [  2.71000000e+04],
        [  5.82330000e+04],
        [  1.22618000e+05],
        [  1.70403000e+05]],

       ・・・

       [[  2.01700000e+03],
        [  4.80000000e+01],
        [  2.40700000e+03],
        [  2.58800000e+03],
        [  3.79900000e+03],
        [  7.28000000e+03]],

       [[  2.01700000e+03],
        [  4.90000000e+01],
        [  2.58800000e+03],
        [  3.79900000e+03],
        [  7.28000000e+03],
        [  1.27850000e+04]]])
ラベルデータ（labels 変数の内容）例
array([[  1.70403000e+05],
       [  1.51829000e+05],
       ・・・
       [  1.27850000e+04],
       [  2.01270000e+04]])

学習と予測処理

発症報告数と週の値では桁数が大きく異なるケースがあるので、このままで効果的な学習ができるかは微妙です。

そこで、ここではバッチ正規化（BatchNormalization）の層を先頭に追加する事で対処してみました。

バッチ正規化によりミニバッチ単位でデータの正規化が行われるようになります。 この場合、ミニバッチのサイズが重要になってくると思います。

また、IDWR では ISO 週番号を使っているようなので、年によって 52週の場合と 53週の場合があります。

ここでは isocalendar でその年の最終週を取得する事で対応しました。（12/28 を含む週がその年の最終週）

sample2.py
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from datetime import datetime
from functools import reduce
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.layers.normalization import BatchNormalization
from keras.layers.recurrent import GRU
from keras.optimizers import Nadam

data_file = sys.argv[1]
item_name = sys.argv[2]
dest_file = sys.argv[3]
predict_size = int(sys.argv[4])
batch_size = int(sys.argv[5])

t = 4
epoch = 5000
n_num = 80

df = pd.read_csv(data_file, encoding = 'UTF-8')

ds = df.groupby(['year', 'week'])[item_name].sum()

def window_with_index(d, wsize):
    return [
        np.r_[
            d.index[i + wsize - 1][0], # 年
            d.index[i + wsize - 1][1], # 週
            d[i:(i + wsize)].values.flatten()
        ] for i in range(len(d) - wsize + 1)
    ]

dw = window_with_index(ds.astype('float'), t + 1)

# 入力データの要素数
input_num = len(dw[0]) - 1

data = np.array([i[0:-1] for i in dw]).reshape(len(dw), input_num, 1)
labels = np.array([i[-1] for i in dw]).reshape(len(dw), 1)

model = Sequential()

# バッチ正規化
model.add(BatchNormalization(axis = 1, input_shape = (input_num, 1)))

model.add(GRU(n_num, activation = 'relu'))

model.add(Dense(1))
model.add(Activation('linear'))

opt = Nadam()

model.compile(loss = 'mean_squared_error', optimizer = opt)

hist = model.fit(data, labels, epochs = epoch, batch_size = batch_size)

fig, axes = plt.subplots(1, 2, figsize = (12, 6))

axes[0].plot(hist.history['loss'])

# 実データの描画
axes[1].plot(ds.values, label = 'actual')

# predict1（入力データをそのまま使用）
res1 = model.predict(data)

axes[1].plot(range(t, len(res1) + t), res1, label = 'predict1')

# 指定年の最終週（52 か 53）の取得
def weeks_of_year(year):
    return datetime(year, 12, 28).isocalendar()[1]

# predict2 用の予測処理
def predict(a, b):
    r = model.predict(a[1])

    year = a[1][0, 0, 0]
    week = a[1][0, 1, 0] + 1

    if week > weeks_of_year(int(year)):
        year += 1
        week = 1

    next = np.r_[
        year,
        week,
        a[1][:, 3:].flatten(),
        r.flatten()
    ].reshape(a[1].shape)

    return (np.append(a[0], r), next)

fst_data = data[0].reshape(1, input_num, 1)
# predict2（入力データの先頭のみ使用）
res2, _ = reduce(predict, range(predict_size), (np.array([]), fst_data))

axes[1].plot(range(t, predict_size + t), res2, label = 'predict2')

min_year = min(ds.index.levels[0])
years = range(min_year, min_year + int(predict_size / 52) + 1)

# X軸のラベル設定
axes[1].set_xticklabels(years)
# X軸の目盛設定
axes[1].set_xticks(
    reduce(lambda a, b: a + [a[-1] + weeks_of_year(b)], years[:-1], [0])
)

# 凡例をグラフ外（下部）へ横並びで出力
axes[1].legend(bbox_to_anchor = (1, -0.1), ncol = 3)

fig.subplots_adjust(bottom = 0.15)

plt.savefig(dest_file)

実行

インフルエンザに対して 2014年 5週目から 250週分をバッチサイズ 52 で予測（predict2）するように実行してみます。

実行例（インフルエンザ）
> python sample2.py idwr.csv インフルエンザ sample2.png 250 52

・・・
Epoch 4998/5000
202/202 [==============================] - 0s 155us/step - loss: 133048118.4158
Epoch 4999/5000
202/202 [==============================] - 0s 155us/step - loss: 154247298.2970
Epoch 5000/5000
202/202 [==============================] - 0s 232us/step - loss: 171130778.4554
1. インフルエンザ （250週予測, batch_size = 52）

f:id:fits:20180121173637p:plain

predict2 は劇的に改善され、将来の予測値もそれっぽいものになりました。

他の感染症に対しても同じように処理してみると以下のようになりました。

2. 感染性胃腸炎 （250週予測, batch_size = 52）

f:id:fits:20180121173656p:plain

3. 手足口病1 （250週予測, batch_size = 52）

f:id:fits:20180121173714p:plain

2年周期を学習できていないような結果となりました。

そこで、バッチサイズを 104 にして 320週分を予測してみたところ以下のようになりました。

4. 手足口病2 （320週予測, batch_size = 104）

f:id:fits:20180121173731p:plain

5. 水痘1 （250週予測, batch_size = 52）

f:id:fits:20180121173748p:plain

predict2 の差異が目立っています。

そこで、バッチサイズを 26 にしてみたところ改善が見られました。 ちなみに、バッチサイズを大きくしても特に改善は見られませんでした。

6. 水痘2 （250週予測, batch_size = 26）

f:id:fits:20180121173805p:plain

7. 流行性耳下腺炎 （250週予測, batch_size = 52）

f:id:fits:20180121173828p:plain