Technology Topics by Brains

ブレインズテクノロジーの研究開発機関「未来工場」で働くエンジニアが、先端オープン技術、機械学習×データ分析(異常検知、予兆検知)に関する取組みをご紹介します。

Pythonによる時系列データの異常検知

インターン生の松井(B4)です.時系列データの異常検知手法をまとめました.

入門 機械学習による異常検知という本の7章が時系列データの異常検知を扱っています.(本書の内容をまとめたWeb記事もあります.)
www.coronasha.co.jp
この本のサンプルコードはすべてRで書かれているため,Python (+numpy, scikit-learn) で書き直してみました.

後半では,深層学習を用いた時系列データの異常検知手法について,知られている所をまとめました.

k近傍法による異常部位検出

時系列データの異常検知手法の中でも比較的シンプルなやり方です.
訓練用の時系列データをスライド窓によりベクトル化しておき,k近傍法を用いて新たな時系列の異常度を計算します.

k=1として最近傍法を実装したコードがこちら.

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import NearestNeighbors

def main():
    data = np.loadtxt("your path/qtdbsel102.txt",delimiter="\t")

    train_data = data[1:3000, 2]
    test_data = data[3001:6000, 2]

    width = 100
    nk = 1

    train = embed(train_data, width)
    test = embed(test_data, width)

    neigh = NearestNeighbors(n_neighbors=nk)
    neigh.fit(train)
    d = neigh.kneighbors(test)[0]

    # 距離をmax1にするデータ整形
    mx = np.max(d)
    d = d / mx

    # プロット
    test_for_plot = data[3001+width:6000, 2]
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax2 = ax1.twinx()

    p1, = ax1.plot(d, '-b')
    ax1.set_ylabel('distance')
    ax1.set_ylim(0, 1.2)
    p2, = ax2.plot(test_for_plot, '-g')
    ax2.set_ylabel('original')
    ax2.set_ylim(0, 12.0)
    plt.title("Nearest Neighbors")
    ax1.legend([p1, p2], ["distance", "original"])
    plt.savefig('./results/knn.png')
    plt.show()


def embed(lst, dim):
    emb = np.empty((0,dim), float)
    for i in range(lst.size - dim + 1):
        tmp = np.array(lst[i:i+dim])[::-1].reshape((1,-1)) 
        emb = np.append( emb, tmp, axis=0)
    return emb

if __name__ == '__main__':
    main()

心電図のデータをサンプルとして用いています.
周期的で正常な心電が見られる部分を訓練データとし,乱れている部分を含むようにテストデータの区間を取りました.
結果のプロットは以下のようになりました.
f:id:Nori_matsu:20170930222409p:plain
異常を含む区間で異常度 (グラフ中distanceと表記) が高くなっていることがわかります.

k近傍法では,注目波形パターンが,訓練データ中の波形パターン群の中でk番目に近いものと,どれだけ異なるかを表しています.
この手法は,綺麗な周期的データが予測される場合には有用ですが,周期的に繰り返しつつも長期的に上昇していくようなトレンドを持つ波形の場合にうまくいきません.

特異スペクトル変換法による変化点検知

一方,ある点の前後の波形を比較して変化点検出を行う特異スペクトル変換法という手法は,波形の長期的なトレンドにも強いと考えられます.
特異スペクトル変換法については,こちらのWebサイトにわかりやすくまとまっています.
numpyを用いて特異スペクトル変換法を実装したコードはこちら.

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

def main():
    data = np.loadtxt("your path/qtdbsel102.txt",delimiter="\t")

    train_data = data[1:3000, 2]
    test_data = data[3001:6000, 2]

    w = 50 # width
    m = 2
    k = w/2
    L = k/2 # lag
    Tt = test_data.size
    score = np.zeros(Tt)

    for t in range(w+k, Tt-L+1+1):
        tstart = t-w-k+1
        tend = t-1
        X1 = embed(test_data[tstart:tend], w).T[::-1, :] # trajectory matrix
        X2 = embed(test_data[(tstart+L):(tend+L)], w).T[::-1, :] # test matrix

        U1, s1, V1 = np.linalg.svd(X1, full_matrices=True)
        U1 = U1[:,0:m]
        U2, s2, V2 = np.linalg.svd(X2, full_matrices=True)
        U2 = U2[:,0:m]

        U, s, V = np.linalg.svd(U1.T.dot(U2), full_matrices=True)
        sig1 = s[0]
        score[t] = 1 - np.square(sig1)

    # 変化度をmax1にするデータ整形
    mx = np.max(score)
    score = score / mx

    # プロット
    test_for_plot = data[3001:6000, 2]
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax2 = ax1.twinx()

    p1, = ax1.plot(score, '-b')
    ax1.set_ylabel('degree of change')
    ax1.set_ylim(0, 1.2)
    p2, = ax2.plot(test_for_plot, '-g')
    ax2.set_ylabel('original')
    ax2.set_ylim(0, 12.0)
    plt.title("Singular Spectrum Transformation")
    ax1.legend([p1, p2], ["degree of change", "original"])
    plt.savefig('./results/sst.png')
    plt.show()


def embed(lst, dim):
    emb = np.empty((0,dim), float)
    for i in range(lst.size - dim + 1):
        tmp = np.array(lst[i:i+dim]).reshape((1,-1))
        emb = np.append( emb, tmp, axis=0)
    return emb

if __name__ == '__main__':
    main()

結果のプロットは以下のようになりました.
f:id:Nori_matsu:20170930230428p:plain
波形が乱れている点において,より鋭く異常度のピークが出ています.

どちらの手法でよりうまく異常検知ができるかはデータの性質に依存しますが,k近傍法の方が計算量は少なくて済むようです.

深層学習を用いた異常検知手法

LSTM (Long short-term memory) を用いた手法

時系列データを深層学習させる手法としては,RNNや,それを発展させたLSTMが知られています.
RNNで来月の航空会社の乗客数を予測する
わかるLSTM ~ 最近の動向と共に

LSTMは,与えられた波形から,次の時刻での値を予測するモデルです.
LSTMを用いて時系列データの異常検知をするアプローチとしては,異常例のデータセットが十分に集められるかそうでないかによって,二種類のアプローチがあります.
異常例のデータセットが少ない場合,直前までの波形から予測された値が,実際の値とどの程度異なるかによって,異常判定ができます.
一方,異常例のデータセットが十分にある場合には,LSTMを直接分類器として用いる方法も用いられるようです.
How to Use LSTM Networks for Time-series Anomaly Detection - Quora

Autoencoder (自己符号化器) を用いた手法

Autoencoderは,ニューラルネットワークを用いて次元圧縮をする手法であり,入力データをそのまま出力 (=復元) させるように教師あり学習をすることで,信号の特徴を抽出することができます.
deepage.net
正常な波形をデータセットとしてAutoencoderを学習させれば,テストデータとして正常波形を用いた時,元の形に近い波形を出力することができます.
一方このモデルに異常な波形を入力した場合,異常な波形に対する特徴を学習していないために,出力波形は大きく異なるでしょう.
このように,部分時系列データをうまく復元できたかどうかによって,異常判定を行います.

結び

時系列データから異常検知を行うためのアルゴリズムとして,k近傍法,特異スペクトル変換法を実装しました.
また,深層学習を用いたアプローチも紹介しました.
Autoencoderによる異常検知アルゴリズムの実装も行いたいところでしたが,時間切れとなってしまったので今回はここまでで...orz
リアルタイムで異常検知を行うためには,精度のみならず,どれだけの計算時間が許されるかという所にも着目して,アルゴリズムを選択する必要があるようです.