はじめに

深層学習を用いた異常検知手法では有名なオートエンコーダを使ってプログラミングしたことをまとめます。オートエンコーダによる再構成誤差とLSTMによる予測誤差などとも比較予定です。

追記
LSTMの異常検知記事、書きました

対象データ

今回もこちらの心電図のデータを使わせていただきました。
心電図のデータですので比較的周期性のあるデータですね。では実装していきます。

なお、データに関しては以前異常検知に関するブログでまとめてあるのでご参照ください。 noleff.hatenablog.com

プログラム

環境

OS: macOS Catalina
numpy: 1.18.5
pandas: 1.0.5
matplotlib: 3.3.0
scikit-learn: 0.23.1
tensorflow: 2.3.0
Keras: 2.4.3

関数

関数定義します。意味はコメントアウトを見てください。

# センサデータ取得
def read_df(path):
    df = pd.read_csv(path)
    return df

# Timeから時間の差を追加
def differ_time(df):
    # Timeから時間の差を割り当てる
    df['dif_sec'] = df['time'].diff().fillna(0)
    df['cum_sec'] = df['dif_sec'].cumsum()
    return df

# 正規化
def act_minxmax_scaler(df):
    scaler = MinMaxScaler()
    scaler.fit(df)
    # 正規化したデータを新規のデータフレームに
    # df_mc = pd.DataFrame(scaler.transform(df), columns=df.columns) 
    # 正規化したデータをnumpyリストに
    mc_list = scaler.fit_transform(df)
    return mc_list

# 部分時系列の作成
def split_part_recurrent_data(data_list, window):
    data_vec = []
    for i in range(len(data_list)-window+1):
        data_vec.append(data_list[i:i+window])
    
    return data_vec

# オートエンコーダの層を作成
def create_auto_encorder(in_out_shape):
    model = Sequential()
    # エンコード
    model.add(Dense(units=200, activation='relu', input_shape=(in_out_shape,)))
    model.add(Dense(units=100, activation='relu'))
    model.add(Dense(units=50, activation='relu'))
    # デコード
    model.add(Dense(units=100, activation='relu'))
    model.add(Dense(units=200, activation='relu'))
    # 出力層
    model.add(Dense(in_out_shape, activation='sigmoid'))

    # 作成したネットワークの確認
    model.summary()

    return model

# オートエンコーダ実行
def act_auto_encoder(model, train_vec, batch_size, epochs):
    # 学習条件の設定　誤差関数＝平均二乗誤差、最適化手法＝Adam法
    model.compile(loss='mse', optimizer='adam', metrics=['acc'])

    # 学習の実行
    hist = model.fit(x=train_vec, y=train_vec, batch_size=batch_size, verbose=1, epochs=epochs, validation_split=0.2)

    return hist

# 評価系のグラフをプロット
def plot_evaluation(eval_dict, col1, col2, xlabel, ylabel):
    plt.plot(eval_dict[col1], label=col1)
    plt.plot(eval_dict[col2], label=col2)
    plt.ylabel(ylabel)
    plt.xlabel(xlabel)
    plt.legend()
    plt.show()

以下ベタ書き

前処理

まずデータを読み込んで、学習データとテストデータにわけます。

# データ読み込み
df = read_df('./ecg.csv')
# 使うデータだけ抽出
test_df = df.iloc[:5000, :] 
train_df = df.iloc[5000:, :]

グラフにするとこんな感じ。

plt.figure(figsize=(16,4))
plt.plot(train_df['signal2'], label='Train')
plt.plot(test_df['signal2'], label='Test')
plt.legend()
plt.show()

f:id:Noleff:20201019182827p:plain

秒数を計算して、時系列情報を秒として持ちます。差分をとったdif_secを加算して、cum_secカラムを新規追加します。

# 秒数を計算
test_df = differ_time(test_df)  
train_df = differ_time(train_df) 

display(test_df)
display(train_df)

f:id:Noleff:20201019182949p:plain

教師なしで異常検知する上で最終的に評価したいので正解のラベルを振ります。
今回は完全にグラフを目で見て、勝手に決めました。
そのため、これが実際に異常かどうかはわかりません。注意してください。
テストデータのグラフを15秒から20秒の間でクローズアップしました。
赤色に示すとおり、17秒から17.5秒の間を異常としました。

anomaly_df = test_df[(test_df['cum_sec']>15)&(test_df['cum_sec']<20)]
plt.figure(figsize=(16,4))
plt.plot(anomaly_df['cum_sec'], anomaly_df['signal2'], color='tab:orange', label='Test')

anomaly_df = test_df[(test_df['cum_sec']>17)&(test_df['cum_sec']<17.5)] 
plt.plot(anomaly_df['cum_sec'], anomaly_df['signal2'], color='tab:red', label='Test')
plt.xlabel('sec')

plt.show()

f:id:Noleff:20201019183505p:plain

赤色の部分にラベル1、それ以外をラベル0にして、新規カラムlabelを定義します。

test_df.loc[(test_df['cum_sec']>17)&(test_df['cum_sec']<17.5), 'label'] = 1
test_df['label'] = test_df['label'].fillna(0)

次に部分時系列を作ってきます。異常検知には「外れ値検知」「変化点検知」「異常部位検知」の大きく3種類ががあります。 3種類の異常に関してはこちらを参照ください。

www.albert2005.co.jp

今回の心電図データは「異常部位検知」に該当します。
異常部位検知において、言葉の通り、ある時系列の部位を異常として検知したいです。そこで時系列を部分ごとに分割した部分時系列を作ります。
生データそのままだとうまく学習できないため、正規化してデータを0~1の範囲に値を収めます。

window = 250 # 部分窓

# 正規化
train_vec = act_minxmax_scaler(train_df[['signal2']])
test_vec = act_minxmax_scaler(test_df[['signal2']])

# 一次元リスト
train_vec = train_vec.flatten()
test_vec = test_vec.flatten()

# 部分時系列作成
train_vec = split_part_recurrent_data(train_vec, window)
test_vec = split_part_recurrent_data(test_vec, window)

学習

続いてメインのオートエンコーダです。学習するための層を構築していきます。
in_out_shapeはオートエンコーダの入力と出力のデータ数です。今回の場合はwindowと同じです。
オートエンコーダの中間層は適当に決めています。ここをチューニングすることこそが、ニューラルネットワークの醍醐味といえるでしょう。

in_out_shape = window
# オートエンコーダ作成
model = create_auto_encorder(in_out_shape)

Model: "sequential"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
dense (Dense)                (None, 200)               50200     
_________________________________________________________________
dense_1 (Dense)              (None, 100)               20100     
_________________________________________________________________
dense_2 (Dense)              (None, 50)                5050      
_________________________________________________________________
dense_3 (Dense)              (None, 100)               5100      
_________________________________________________________________
dense_4 (Dense)              (None, 200)               20200     
_________________________________________________________________
dense_5 (Dense)              (None, 250)               50250     
=================================================================
Total params: 150,900
Trainable params: 150,900
Non-trainable params: 0
_________________________________________________________________

学習を実行させていきます。オートエンコーダは学習データとして入力したデータを再現するように学習するニューラルネットワークです。
そのため、説明変数と目的変数は同じtrain_vecが入ります（act_auto_encoder関数の話）。

batch_size = 100
epochs = 20

train_vec = np.array(train_vec)
test_vec = np.array(test_vec)
hist = act_auto_encoder(model, train_vec, batch_size, epochs)

Epoch 1/20
318/318 [==============================] - 1s 3ms/step - loss: 0.0090 - acc: 0.0584 - val_loss: 0.0014 - val_acc: 0.1011
Epoch 2/20
318/318 [==============================] - 1s 3ms/step - loss: 0.0011 - acc: 0.1096 - val_loss: 7.4762e-04 - val_acc: 0.0962
Epoch 3/20
318/318 [==============================] - 1s 3ms/step - loss: 6.4285e-04 - acc: 0.1120 - val_loss: 5.4386e-04 - val_acc: 0.0996
Epoch 4/20
318/318 [==============================] - 1s 3ms/step - loss: 5.2786e-04 - acc: 0.1153 - val_loss: 4.9082e-04 - val_acc: 0.1028
Epoch 5/20
318/318 [==============================] - 1s 3ms/step - loss: 4.7851e-04 - acc: 0.1195 - val_loss: 4.2493e-04 - val_acc: 0.1082
Epoch 6/20
318/318 [==============================] - 1s 3ms/step - loss: 4.3090e-04 - acc: 0.1236 - val_loss: 3.7944e-04 - val_acc: 0.1099
Epoch 7/20
318/318 [==============================] - 1s 3ms/step - loss: 4.0277e-04 - acc: 0.1275 - val_loss: 3.7550e-04 - val_acc: 0.1146
Epoch 8/20
318/318 [==============================] - 1s 3ms/step - loss: 3.8248e-04 - acc: 0.1302 - val_loss: 3.4081e-04 - val_acc: 0.1204
Epoch 9/20
318/318 [==============================] - 1s 3ms/step - loss: 3.7038e-04 - acc: 0.1336 - val_loss: 3.3807e-04 - val_acc: 0.1229
Epoch 10/20
318/318 [==============================] - 1s 3ms/step - loss: 3.5362e-04 - acc: 0.1362 - val_loss: 3.4684e-04 - val_acc: 0.1269
Epoch 11/20
318/318 [==============================] - 1s 3ms/step - loss: 3.4186e-04 - acc: 0.1391 - val_loss: 3.1406e-04 - val_acc: 0.1294
Epoch 12/20
318/318 [==============================] - 1s 3ms/step - loss: 3.2798e-04 - acc: 0.1446 - val_loss: 3.1423e-04 - val_acc: 0.1332
Epoch 13/20
318/318 [==============================] - 1s 3ms/step - loss: 3.2826e-04 - acc: 0.1473 - val_loss: 2.9715e-04 - val_acc: 0.1377
Epoch 14/20
318/318 [==============================] - 1s 3ms/step - loss: 3.1464e-04 - acc: 0.1475 - val_loss: 3.0181e-04 - val_acc: 0.1439
Epoch 15/20
318/318 [==============================] - 1s 3ms/step - loss: 3.0473e-04 - acc: 0.1517 - val_loss: 2.9136e-04 - val_acc: 0.1446
Epoch 16/20
318/318 [==============================] - 1s 4ms/step - loss: 2.9575e-04 - acc: 0.1564 - val_loss: 2.9008e-04 - val_acc: 0.1465
Epoch 17/20
318/318 [==============================] - 1s 3ms/step - loss: 2.9276e-04 - acc: 0.1569 - val_loss: 2.8658e-04 - val_acc: 0.1494
Epoch 18/20
318/318 [==============================] - 1s 3ms/step - loss: 2.8440e-04 - acc: 0.1574 - val_loss: 2.7905e-04 - val_acc: 0.1521
Epoch 19/20
318/318 [==============================] - 1s 3ms/step - loss: 2.7920e-04 - acc: 0.1611 - val_loss: 2.6296e-04 - val_acc: 0.1561
Epoch 20/20
318/318 [==============================] - 1s 4ms/step - loss: 2.7432e-04 - acc: 0.1648 - val_loss: 2.5073e-04 - val_acc: 0.1508

戻り値histから、誤差値と評価値の収束具合を描画します。
val_lossよりもlossが大きければ未学習の可能性が高いです。

# 誤差の収束具合を描画
plot_evaluation(hist.history, 'loss', 'val_loss', 'epoch', 'loss')

f:id:Noleff:20201019184351p:plain

val_accよりもaccが大きければ過学習の可能性が高いです。

# 評価の収束具合を描画
plot_evaluation(hist.history, 'acc', 'val_acc', 'epoch', 'acc')

f:id:Noleff:20201019184404p:plain

異常検知

作成したモデルから予測してみます。
異常部位の予測が変になってますね。

# モデルをテストデータに適用
pred_vec = model.predict(test_vec)

# テストデータと出力データの可視化
plt.plot(test_vec[:,0], label='test')
plt.plot(pred_vec[:,0], alpha=0.5, label='pred')
plt.legend()
plt.show()

f:id:Noleff:20201019184456p:plain

テストデータとテストデータを再現したデータの再構成誤差を計算します。シンプルな距離の計算です。
グラフをみると、まあまあうまくできているように見えます。

# 再構成誤差（異常スコア）の計算
dist = test_vec[:,0] - pred_vec[:,0]
dist = pow(dist, 2)
dist = dist / np.max(dist)

# 異常スコアの可視化
plt.plot(dist)
plt.show()

f:id:Noleff:20201019184709p:plain

もとのデータと比べてみるとこんな感じです。

plt.plot(test_vec[3800:,0], label='test')
plt.plot(pred_vec[3800:,0], label='pred')
plt.plot(dist[3800:], label='dist')
plt.legend()
plt.ylabel('anomaly score')
plt.show()

f:id:Noleff:20201019190850p:plain

評価

実際に評価してみます。まず、異常度からラベルが0か1かを決めます。異常度は0から1に収まりますので、試しに0.8以上を異常としてみます。

pred_label_list = [1 if d >= 0.8 else 0 for d in dist]

pred_df = pd.DataFrame(pred_label_list, columns=['pred_label'])

new_test_df = pd.concat([test_df, pred_df], axis=1)
new_test_df = new_test_df.fillna(0)

0と1の離散型のデータをグラフにしているため、カクカクしてます。
青色が初めに定義した異常のラベル。オレンジ色が異常度が0.8以上を異常としたラベルです。

plt.plot(new_test_df['label'], label='true anomaly')
plt.plot(new_test_df['pred_label'], label='pred anomaly')
plt.ylabel('anomary label')
plt.legend()

f:id:Noleff:20201019193921p:plain

分類タスクで評価に使われる混同行列から、精度、適合率、再現率を計算します。

print(confusion_matrix(new_test_df['label'], new_test_df['pred_label']))
print('accuracy_score: ', accuracy_score(new_test_df['label'], new_test_df['pred_label']))
print('precision_score: ', precision_score(new_test_df['label'],new_test_df['pred_label']))
print('recall_score: ', recall_score(new_test_df['label'],new_test_df['pred_label']))

果たしてこれはうまくいっていると言えるでしょうか。
結果から分かる通り、異常の誤検知は一回もしていません。それは適合率から明らかです。
ですが、再現率が非常に低いことから、多くの異常を見逃していることがわかります。

[[4876    0]
 [ 110   14]]
accuracy_score:  0.978
precision_score:  1.0
recall_score:  0.11290322580645161

ここで異常度からラベルを振る閾値をまじめに考えてみます。
まずは異常度をヒストグラムで描画してみます。
ほとんど0近辺の値しかないことがわかります。
つまり、閾値0.8は本当にかけ離れた異常度のみを異常としてラベルをふっていることになります。

f:id:Noleff:20201019190928p:plain

閾値を0.1 以上にしてみましょう。
閾値を下げたのでもちろん誤検知の数は増えます。しかし、同時に異常の見逃しの数を減らすことができました。
誤検知と見逃しはトレードオフの関係なので、どういう異常検知をしたいかで閾値を調整する必要があります。

pred_label_list = [1 if d >= 0.1 else 0 for d in dist]

pred_df = pd.DataFrame(pred_label_list, columns=['pred_label'])

new_test_df = pd.concat([test_df, pred_df], axis=1)
new_test_df = new_test_df.fillna(0)

plt.plot(new_test_df['label'], label='true anomaly')
plt.plot(new_test_df['pred_label'], label='pred anomaly')
plt.ylabel('anomary label')
plt.legend()
plt.savefig('figure/11_0.1以上の異常ラベル.png')

print(confusion_matrix(new_test_df['label'], new_test_df['pred_label']))
print('accuracy_score: ', accuracy_score(new_test_df['label'], new_test_df['pred_label']))
print('precision_score: ', precision_score(new_test_df['label'],new_test_df['pred_label']))
print('recall_score: ', recall_score(new_test_df['label'],new_test_df['pred_label']))

[[4852   24]
 [  46   78]]
accuracy_score:  0.986
precision_score:  0.7647058823529411
recall_score:  0.6290322580645161

f:id:Noleff:20201019191710p:plain

まとめ

オートエンコーダによる時系列データの教師なし異常検知をプログラミングしました。
時系列データは有名な心電図のデータです。
簡単にプログラムしたため、異常検知の精度はいまいちです。心電図の生データを正規化しただけの前処理しかしていないため、もう少し工夫をする必要があるかもしれません。
また、オートエンコーダの層の次元数に関しても色々試してみる価値はあると思います。