【時系列データ】深層距離学習で異常検知

以前に、[【まとめ】ディープラーニングを使った異常検知](https://qiita.com/shinmura0/items/06d81c72601c7578c6d3)の記事を
書きましたが、主に「画像」を対象とした内容でした。

今回は、その記事で有効だった深層距離学習を「時系列データ」に
適用して異常検知を行います。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264781/40111808-8e8a-8bb3-a414-849cad28f883.png)

#画像データとの違い
画像の場合は、比較的有利な条件で学習が可能でした。
一方、**時系列データの場合は**、以下の制限があります。

+ 転移学習が使えない
画像では学習済モデルが公開されているため、転移学習が使えましたが、
時系列データの学習済モデルは一般的に公開されていないため、転移学習が使えません。

+ 参照データがない
metric learningを使う場合、正常データと共に全然関係ない参照データ
（正常データと見比べるデータ）が必要でした。ところが、時系列データの
場合、参照データが存在しない可能性が高いです。

以上の制約から、時系列データは自作のモデル+正常データのみで学習させる
必要があります。そこで、本稿では**自己教師あり学習**を使います。

#自己教師あり学習
本稿では、正常データに強引にラベルを付けて、自己教師あり学習を行います。
ラベルを付けて深層距離学習を適用することにより、時系列データの特徴を捉えることが
でき、最終的に異常検知が可能となります。

本稿では、ラベルはクラスタリングで付けていきます。

時系列データのクラスタリングは様々な手法が提案されていますが、ここでは
k-means法を使います。

#実装
[カルマンフィルタの記事](https://qiita.com/shinmura0/items/8c1e83636291f10e0840)と同じく、心拍数のデータを使います。

##データロード

```python
import urllib.request
import numpy as np

import os
import keras
import matplotlib.pyplot as plt

from keras.utils import to_categorical
from keras.layers import Input, GlobalAveragePooling2D, Dense, Dropout
from keras.layers import Activation, BatchNormalization, Conv2D
from keras.initializers import he_normal
from keras.optimizers import Adam
from keras.models import Model
from keras import backend as K

from sklearn.preprocessing import MinMaxScaler
from sklearn.cluster import KMeans
from sklearn.neighbors import LocalOutlierFactor

url =  "http://www.cs.ucr.edu/~eamonn/discords/qtdbsel102.txt"
f_name = "data.txt"
 
# ダウンロードを実行
urllib.request.urlretrieve(url, f_name)

# データ抽出
data = np.loadtxt(f_name, delimiter="\t")
data = data[2000:5000,1]

plt.plot(data)
plt.show()
```

0～1000は学習データ、1000～3000はテストデータで使用します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264781/67260d91-3ae7-0742-c517-357643684601.png)


##時系列データを作成
検出窓を一つずつズラしながらデータを作成します。

```python
w = 50 # 抽出窓の幅
T = 1000 # 学習データの範囲

x, x_test = [], []
    
for i in range(w, T):
    x.append(data[i-w: i])

for i in range(T, len(data)-w):
    x_test.append(data[i: i+w])

x = np.array(x)
x_test = np.array(x_test)

#0-1
x_test = (x_test - np.min(x))/(np.max(x) - np.min(x))
x = (x - np.min(x))/(np.max(x) - np.min(x))
```

##k-means法
ここでは、クラスターの数を5個にしています。

```python
classes = 5
y_train = KMeans(n_clusters=classes).fit_predict(x)
Y_train = to_categorical(y_train)
print(np.sum(Y_train, axis=0))
```

```
[ 79. 428. 248.  95. 100.]
```

かなり不均衡なデータ数になりましたが、このままCNNに突っ込みます。
（ノイズを入れながら、DataAugmentationして均衡にしてみましたが、
結果はさほど変わりませんでした。）

##深層距離学習
L2-SoftmaxLossを使って学習させます。

```python
def train(x, y):
    #cnnの構築
    alpha = 5

    inputs = Input(shape=x.shape[1:])

    c = Conv2D(64, (1, 1), padding="same", kernel_initializer=he_normal())(inputs)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Conv2D(64, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Dropout(0.2)(c)
    c = Conv2D(64, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Conv2D(64, (1, 3), strides=2, kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)

    c = Conv2D(128, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Conv2D(128, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Dropout(0.2)(c)
    c = Conv2D(128, (1, 3), strides=2, kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)

    c = Conv2D(256, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Conv2D(256, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Dropout(0.2)(c)
    c = Conv2D(256, (1, 3), strides=2, kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)

    c = Conv2D(512, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Conv2D(512, (1, 3), padding="same", kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu")(c)
    c = Dropout(0.2)(c)
    c = Conv2D(512, (1, 3), strides=2, kernel_initializer=he_normal())(c)
    c = BatchNormalization()(c)
    c = Activation("relu", name="out")(c)

    c = GlobalAveragePooling2D()(c)
    c = keras.layers.Lambda(lambda xx: alpha*(xx)/K.sqrt(K.sum(xx**2)))(c) #metric learning
    c = Dense(y.shape[1], activation='softmax')(c)
    model = Model(inputs, c)
    #model.summary()

    model.compile(loss='categorical_crossentropy',
                  optimizer=Adam(lr=0.0001, amsgrad=True),
                  metrics=['accuracy'])

    #cnnの学習
    hist = model.fit(x, y, batch_size=32, epochs=50, verbose=True)

    return model
```

データを4次元にしてから学習させます。

```python
X_train = np.expand_dims(x, axis=1)
X_test = np.expand_dims(x_test, axis=1)

X_train = np.expand_dims(X_train, axis=-1)
X_test = np.expand_dims(X_test, axis=-1)

model = train(X_train, Y_train)
```

##結果
最後にLOFを使ってテストデータの異常値を取得します。

```python
def get_score(model, x_train, x_test):
    model_s = Model(inputs=model.input,outputs=model.layers[-2].output)
    train = model_s.predict(x_train, batch_size=1)
    test = model_s.predict(x_test, batch_size=1)
    
    train = train.reshape((len(x_train),-1))
    test = test.reshape((len(x_test),-1))
    
    ms = MinMaxScaler()
    train = ms.fit_transform(train)
    test = ms.transform(test)
    
    # fit the model
    lof = LocalOutlierFactor(n_neighbors=5)
    y_pred = lof.fit(train)

    # plot the level sets of the decision function
    Z = -lof._decision_function(test)

    return Z
```

```python
score = get_score(model, X_train, X_test)
```

```python
plt.figure(figsize=(12,12))
plt.subplot(2,1,1)
plt.plot(data[1000:], label="true")
plt.legend()

plt.subplot(2,1,2)
plt.plot(np.arange(50,2000),score, label="score")
plt.legend()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264781/a8c1483f-4055-faf5-0d5e-9149dd509b86.png)

ちゃんと異常部分が検知できました。
参考までにt-sneで可視化してみます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264781/40111808-8e8a-8bb3-a414-849cad28f883.png)

黒い点が学習データ、カラフルな点がテストデータです。

まず、今回はクラスター数を5にしており、大体5つのグループに分かれて
いることが分かります。（黒い丸で囲ってみました。）

さらに、テストデータでは、時間1250付近に異常な波形が見られます。
その付近（青い点）を注視してみると、赤い矢印で示した点だけが仲間外れに
なっており、「距離」学習がちゃんと機能していることが分かります。


#考察
+ 再現性が低いため、実際の運用では多数のCNNモデルを用意して、アンサンブルを取った方が良さそう。
+ このモデル自体はRNNやカルマンフィルタなど、時系列で情報の引き継ぎはしていないため、時系列データのみならず、テーブルデータなどでも効果があると思われる。

##まとめ
+ はっきり言って、（カルマンフィルタなどと比べ）計算量が重いため、**あえて使うメリットはない**
+ ただ、使い方を覚えてしまえば、「画像」でも「時系列データ」でも処理してしまう「万能な検知器」として使える
+ 従来手法（K近傍法やカルマンフィルタなど）を使っても**うまくいかない場合**、突破口の一つになるかもしれない