初心者のための時系列データ異常検知

# はじめに
この記事は、シスコシステムズ合同会社の同士による [Cisco Advent Calendar 2019](https://qiita.com/advent-calendar/2019/cisco) の 6 日目として投稿しています。

2012 年に「ImageNet Large Scale Visual Recognition Challenge（ILSVRC） 2012」でカナダ・トロント大学の SuperVision チームが圧勝してから、7 年が経ちました。この時、深層学習を用いたモデルを用いることによって、第 3 次 AI ブームでは機械学習や深層学習でのモデリングが多く行われ、昨今では、AI と言えば機械学習・深層学習を用いることが常識となっています。

しかしながら、2017 年に発表された総務省のアンケートによると、AI ソリューションを実際に導入した企業は、14.1% にとどまっており、導入を検討している企業も 22.8% であるなど、AI ソリューションに前向きな企業は検討段階を含めても 36.9% という状況です。
（引用：[AI・IoT導入状況と予定](http://www.soumu.go.jp/johotsusintokei/whitepaper/ja/h30/html/nd132210.html)　総務省）

機械学習に関するオープンソースのサービスは大変多く、小規模での利用も容易です。そんな中、過去 2 年に渡って、私は機械学習を簡単に利用するための例を記事として書かせていただきました。
[2017 年記事： API を使って得た位置情報を機械学習させてみた](https://qiita.com/kmasuyam/items/8aaf54f08ae87f5646ce)
[2018 年記事： コラボレーションツールを使った機械学習教師データ集め](https://qiita.com/kmasuyam/items/7a1c020c9c600a4eeae3)

そこで今回は、データとして 1 番多く取得されるであろう時系列データを使って、簡単に機械学習のモデリングができる例をご紹介したいと思います。その中でも、異常検知という分野ではどのように機械学習が応用されているかにフォーカスします。

# 時系列データ
時系列データとは、時間の推移とともに観測されるデータ全般のことを指します。そして、観測される順序に意味を持つものです。Cisco 製品に関わる部分での時系列データですと、Syslog データも時系列データに含まれますし、[Stealthwatch](https://www.cisco.com/c/ja_jp/products/security/stealthwatch/index.html) の Flow データや、[Cloudlock](https://www.cisco.com/c/ja_jp/products/security/cloudlock/index.html) の接続先データなど、多岐に渡ります。

時系列データを機械学習・深層学習に応用する場合は、下記のようなフローをとることが多くなります。
![4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214710/a852f07b-de89-5764-14af-0d9b91b5fede.png)
類似度に基づく手法では、特徴量をどのように抽出するか、時系列データのウィンドウサイズ（特徴量が現れやすい時系列セットの区切り方）の設定、後の機械学習手法におけるモデルをどのように選択するかによって、その精度が大きく変化します。

モデルに基づく手法では、まず計量時系列分析に用いられる手法を対応させます。ここで言う計量時系列分析に用いられる手法とは、隠れマルコフモデルや ARMA モデルや、VAR モデル、SARIMA モデルなどのことを指します。状態遷移を時系列の変化と見做し、データの遷移を表すモデルのことです。これらの手法では、各モデルにパラメータを持つため、得られたパラメータを用いて機械学習手法に適用します。


<img width="462" alt="スクリーンショット 2019-12-06 18.31.24.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214710/fea0453f-f8d5-9620-5175-27e9ab636d22.png">
図 1: Cisco Web Security Appliance の時系列データ

<img width="1031" alt="スクリーンショット 2019-12-06 18.32.00.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214710/5f639f5f-d505-11c0-8a57-59429cca49dc.png">
図 2: Cisco Cloudlock の時系列データ


今回は、数値データとして得られた時系列データに対して、「通常のデータと比較し、異常なデータを発見する」ことを目標にします。


# 異常検知
異常検知は、統計学、確率論、最適化論など、数式が並ぶ大変難しい分野だと私は認識していますが、適用できる分野は非常に多く
- 工場での異常兆候発見
- Security 分野における Malware 等の問題発見
- 人体測定データを利用した健康に関する異常発見
など、さまざまな分野が考えられます。これらの分野で異常検知ができることで、人間が問題を検知するのに比べて早急に対処できるなどのメリットが得られます。

異常検知の分野では、多くの機械学習手法がその適用方法を提案されていますが、ここでは一部をご紹介します。


異常検知の分野では、重要視する点が 2 通り存在する、と個人的には認識しています。それは、精度 (特に False-Negative を減らす考え方) を重要視する考え方と、検知までのスピードを重要視する考え方です。これは、実際の運用をどのように行うかにも関わるため、重要な考え方です。

### 1. 異常検知の精度を重要視する考え方
精度を重要視する考え方は、特に故障検知の分野に多く見られます。これは、```本当は故障しているのに、「故障していない」という予測が出る```ことを避けたいためです。
機械学習分野の用語で言えば、再現率 (recall) を上げる考え方です。

（参考） 再現率とは
精度を考える際は、下記のような混合行列を考慮します。
![5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214710/93f705ac-08f9-d33b-ff91-4aaff057d7f2.png)
図 3: 混合行列

そして、``再現率 (recall)`` とは下記の式で表される精度のことです。
![6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214710/6c0f5ce1-4192-6f84-6e8e-159a27dcc78f.png)
つまり、再現率とは「実際には異常であるときに、モデルがどの程度異常と判定できるか」という指標です。

現実的には、工場などの設備が故障しているかを判定するのに、人間が関わらない運用は難しいと考えられます。しかし、False-Negative の数が多い場合、AI アプリケーションへの信用が少ないために導入する前後で運用が変わらず、管理者の負担が減らないことが予想されるのです。
この状況を避けるために、False-Negative の数を減らし、故障検知の負担を減らすような AI アプリケーションが期待されます。

### 2. 検知までのスピードを重要視する考え方
機械学習・深層学習を利用する際は、計算コストも考えなければなりません。モデルの精度がいくら高くても、検知するまでにかかる時間が遅い場合、AI アプリケーションを導入する意味がないからです。

例えば、セキュリティに関する AI アプリケーションがこの場合に当たります。いくら異常なトラフィックを検知したとしても、検知されたのが翌日の場合、すでに重要な情報を盗まれた後で、何の意味も成しません。
つまり、サーバなどの基盤に大金を投入できる場合以外は、モデルの計算コストを考慮し、検知までのスピードを重要視する必要があります。

計算コストが低い機械学習手法には、単純ベイズ分類器や k-means/k-medoids 法などの手法が挙げられます。ここでは詳細な手法の説明は避けますが、この他にも異常検知で用いられる手法は数多く存在します。

# 適用例
ここでは、適用例を考えてみたいと思います。利用するデータは、Stealthwatch の Flow 数のような 2 次元時系列データを想定します。できる限り一般化したいため、データは想定とします。
<img width="1427" alt="スクリーンショット 2019-12-06 18.31.42.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214710/de5ca695-5dc0-98f5-87ea-f1fc0f72e352.png">
図 4: Cisco Stealthwatch の時系列データ

上記のような時系列データから、ウィンドウサイズ分のデータをパターンとして抽出し、通常状態のパターンと比較します。

まず、パターンの抽出ですが、今回考えなければならないのは、Flow 数が 0 になることはほとんどなく、特にウィンドウサイズとなる数分〜数時間の区切りでは常にトラフィックが流れていると考えられることです。
数値が 0 となることが分かっている場合は、数値が 0 となる時点で区切れば良いのですが、数値が 0 とならない場合は、パターンの抽出を考えなければなりません。

ウィンドウサイズが小さい場合、同じパターンが現れる時間帯や時間幅が、状況によって異なることが予想されます。例えば、社員が昼休憩中に YouTube を見ることが分かっているため、12〜13 時に Flow 数が急激に上がり、13 時を過ぎると分かっていても、日によって最大値をとる時間は異なると考えられるからです （12:31 に最大値をとるか、12:36 に最大値をとるか分からない）。

このように、時間幅に関わらず「形が似ている」場合に、その似ている度合いを算出するのが、Dynamic Time Warping - DTW という手法です。この手法では、データ長が揃っていなくても判定することができます。つまり、時間幅が異なる場合でも判定することができるのです。

反対に、時間軸に対してのズレが重要な意味を持つと考えられる場合、ユークリッド距離を用います。例えば、ウィンドウサイズが大きい場合における、昼間と夜間の Flow 数の違いなどです。

それでは、実際に DTW を Python で実装してみます。DTW の実装には、動的計画法を用います。
また、2 点間の距離を測定する場合に、絶対値で計測する場合とユークリッド距離で計測する場合があると考えられるため、``method`` の引数で分けています。
第 1 引数と第 2 引数は、ウィンドウサイズごとに区切った時系列データ 2 つを表します。数日分の時系列データで比較する場合は、この関数を複数回呼び出す処理をすることになるわけです。

```Python
def dtw(wave_x, wave_y, method="abs"):
    d = np.zeros([len(wave_x)+1, len(wave_y)+1])
    d[:] = np.inf
    d[0, 0] = 0
    if method = "euclid":
        for i in range(1, d.shape[0]):
            for j in range(1, d.shape[1]):
                cost = np.sqrt((wave_x[i-1] - wave_y[j-1])**2)
                cost = (wave_x[i-1] - wave_y[j-1])
                row.append(cost)
                d[i, j] = cost + min(d[i-1, j], d[i, j-1], d[i-1, j-1])
    else:
        for i in range(1, d.shape[0]):
            for j in range(1, d.shape[1]):
                cost = np.abs(wave_x[i-1] - wave_y[j-1])
                row.append(cost)
                d[i, j] = cost + min(d[i-1, j], d[i, j-1], d[i-1, j-1])
    elapsed_time = time.time() - start_time
    return d[-1][-1], d, matrix
```

この DTW を複数回計算することで、複数の時系列データ間の距離行列が求められます。この距離行列を用いて、k-medoids 法による分類を考えます。今回は、「通常/異常」の 2 パターンの分類を目標としているため、クラスタ数を 2 と設定しています。

```Python
self.n_cluster = 2
```

k-medoids 法の実装は下記です。scikit-learn には k-medoids 法が実装されていないため、下記のように実装します。``D_matrix`` の部分には、距離行列を代入します。

```Python
class KMedoids():
    def __init__(self, max_iter=300):
        self.n_cluster = 2
        self.max_iter = max_iter

    def fit_predict(self, D_matrix):
        m, n = D_matrix.shape
        ini_medoids = np.random.choice(range(m), self.n_cluster, replace=False)
        tmp_D = D_matrix[:, ini_medoids]

        labels = np.argmin(tmp_D, axis=1)

        results = pd.DataFrame([range(m), labels]).T
        results.columns = ['id', 'label']

        col_names = ['x_' + str(i + 1) for i in range(m)]
        results = pd.concat([results, pd.DataFrame(D_matrix, columns=col_names)], axis=1)

        old_medoids = ini_medoids
        new_medoids = []

        loop = 0
        while ((len(set(old_medoids).intersection(set(new_medoids))) != self.n_cluster) 
               and (loop < self.max_iter) ):
        if loop > 0:
            old_medoids = new_medoids.copy()
            new_medoids = []
        for i in range(self.n_cluster):
            tmp = results[results['label'] == i].copy()
            tmp['distance'] = np.sum(tmp.loc[:, ['x_' + str(id + 1) for id in tmp['id']]].values, axis=1)
            tmp = tmp.reset_index(drop=True)
            new_medoids.append(tmp.loc[tmp['distance'].idxmin(), 'id'])

        new_medoids = sorted(new_medoids)
        tmp_D = D_matrix[:, new_medoids]

        clustaling_labels = np.argmin(tmp_D, axis=1)
        results['label'] = clustaling_labels
        loop += 1
        results = results.loc[:, ['id', 'label']]
        results['flag_medoid'] = 0

        for medoid in new_medoids:
            results.loc[results['id'] == medoid, 'flag_medoid'] = 1
        tmp_D = pd.DataFrame(tmp_D, columns=['medoid_distance'+str(i) for i in range(self.n_cluster)])
        results = pd.concat([results, tmp_D], axis=1)

        self.results = results
        self.cluster_centers_ = new_medoids
        return results['label'].values
```

上記で 2 クラス分類を行い、通常と異常の 2 つに分類することが可能です。
k-medoids 法の詳細はここでは省きますが、特徴としては k-means 法とあまり変わらないものの、k-means 法と違い medoid を計算して分類するため、外れ値に強いという特徴があります。また、距離行列さえ求められれば分類することが可能なため、応用が効きます。
medoid の求め方は下記です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214710/45ef8e15-fe32-7a6f-b2e4-a697ecfa17a2.png)


# 最後に
今回は、異常検知の手法や適用方法について書かせていただきました。また、その後 2 次元時系列データを用いて異常検知を行なうためのコード例を示しました。

k-medoids 法は、距離行列さえ求められれば分類させることが可能なため、文字列であっても距離を求めることで分類が可能です。文字列の距離は、ジャロ・ウィンクラー距離やレーベンシュタイン距離を用いるため、ぜひこちらも検索してみてください。

ぜひこの記事を参考に、手元にあるデータを使った異常検知を試していただけると嬉しいです。

# 免責事項
本サイトおよび対応するコメントにおいて表明される意見は、投稿者本人の個人的意見であり、シスコの意見ではありません。本サイトの内容は、情報の提供のみを目的として掲載されており、シスコや他の関係者による推奨や表明を目的としたものではありません。各利用者は、本 Web サイトへの掲載により、投稿、リンクその他の方法でアップロードした全ての情報の内容に対して全責任を負い、本 Web サイトの利用に関するあらゆる責任からシスコを免責することに同意したものとします。