テストエンジニアが機械学習してみた備忘録

広島とジト目が好きなテストエンジニアが機械学習に手を出した備忘録。

「入門 機械学習による異常検知」の時系列データの異常検知をPythonでやってみる 〜特異スペクトル変換法編〜

gratk.hatenablog.jp

前回の続きで特異スペクトル変換法による異常検知を試してみます。 特異スペクトル変換法では、特異値分解を行うようなのでその辺りと合わせて整理してみます。

実行環境

  • Mac OS High Sierra (10.13.1)
  • Python 3.6.2
    • numpy 1.13.1
    • matplotlib 2.0.2
    • pandas 0.20.3
    • sklearn 0.19.1
  • jupyter notebook
    • 掲載コードはjupyter notebook上で動作させています

sklearn.decomposition

sklearn.decompositionモジュールは行列分解(matrix decomposition)のためのアルゴリズムを搭載しているとのこと。 主成分分析(principal component analysis; PCA)のためのメソッドも様々なアルゴリズムのものが用意されているようですが、 特異値分解(singular value decomposition; SVD)のためのメソッドも用意されています。

Pythonでの特異値分解

Numpyにも特異値分解のためのメソッドとしてnumpy.linalg.svdが用意されているようです。 scikit-learnとNumpyのどちらを使うかについては、結果はどちらでも変わらないとは思いますが、 以下の記事によるとscikit-learnの方が高速に動作するようです。

soralab.space-ichikawa.com

となるとscikit-learnを用いたいのですが、ちょっと私の前提知識と読解力が足りていないので、 今回は直感的にわかりやすかったNumpyを用います。

numpy.linalg.svd

numpy.linalg.svd — NumPy v1.13 Manual

U, s, V = np.linalg.svd(X, full_matrices=False)

とすることでXを特異値分解し、Uに左特異ベクトルの行列、Vに右特異ベクトルの行列、sに特異値の対角行列を得られます。full_martices=Trueで実行すると、UはXの行数の正方行列、VはXの列数の正方行列となるようです。

やってみる

では事前調査が終わったので実際にやってみます。前回と同じデータを用いますが、訓練データは不要なので、検知したい異常箇所を持つ3001件目〜6000件目のデータを用います。

%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt


'''
dataをsize毎のスライス窓に分割
'''
def embed(data, size):
    window = np.empty((0, size))
    for index in range(0, len(data)-size+1):
        new_window = [data[i] for i in range(index, index+size)]
        window = np.append(window, np.array([new_window]), axis=0)
    return window


df = pd.read_csv("./qtdbsel102.csv", sep="\t", header=None)

# 前から3001件目〜6000件のデータを対象とする
data = df.loc[3000:5999, 2].reset_index(drop=True)

# 窓幅、履歴行列の列サイズ、ラグ、パターン数(書籍中のRでの実行例に従う)
w = 50
K = 25
L = 12
m = 2
Tt = len(data)

abnormality = [0 for i in range(0, Tt)]
for t in range(w+K-1, Tt-L):
    begin_at = t-w-K+1
    end_at = t
    X1 = pd.DataFrame(embed(data[begin_at:end_at].as_matrix(), w)).T
    X1 = X1.iloc[::-1]
    
    begin_at = t-w-K+1+L
    end_at = t+L
    X2 = pd.DataFrame(embed(data[begin_at:end_at].as_matrix(), w)).T
    X2 = X2.iloc[::-1]
    
    # 特異値分解
    U1, s1, V1 = np.linalg.svd(X1, full_matrices=False)
    U1 = U1[:, 0:m]
    U2, s2, V2 = np.linalg.svd(X2, full_matrices=False)
    U2 = U2[:, 0:m]
    
    # 部分空間同士の重なり合いと異常度
    U3, s3, V3 = np.linalg.svd(np.dot(U1.T , U2))
    sig1 = s3[0]
    abnormality[t] = 1 - sig1 * sig1

# テストデータと異常度のプロット
fig, ax1 = plt.subplots()
ax1.plot(test, color='#377eb8')
ax2 = ax1.twinx()
ax2.plot(abnormality, color='#ff7f00')
plt.show()

f:id:gratk:20171214001112p:plain

いい感じにできたんじゃないでしょうか。

まとめ

細かい理屈は全然理解しきれていませんが、ライブラリの力を借りて特異スペクトル変換を試すことはできました。 実務に適用する際は、w, k, L, mといったパラメータを調整するためにもう少し手法に対する理解を深める必要があると思いますが、それは今後の課題としたいと思います。