COMPANY

PAGE TOP

BLOG

RapidMinerとPythonを使用した特異スペクトル変換法による変化検知

皆さん、こんにちは!
KSKアナリティクスの岡地です。

今回は、特異スペクトル変換法を使用した変化検知をRapidMinerとPythonを使用して
実装してみたいと思います。
以下の書籍に特異スペクトル変換法とその実装についての説明があったので参考にしました。
井手剛・杉山将(2015)「異常検知と変化検知(機械学習プロフェッショナルシリーズ)」,講談社.

使用するデータは心電図のデータを使用します。
元データは以下よりダウンロードできます。
http://www.cs.ucr.edu/~eamonn/discords/qtdbsel102.txt
今回は上記よりダウンロードしたデータを加工したデータを使用します。
使用するデータとRapidMinerのプロセスはこちらよりダウンロードして下さい。
※プロセスの実行にはPython Scripting Extensionが必要です。
 メニューバー>Extensions>Marketplaceより無償でインストールできます。

今回のデータは基本的には周期的な動きをしていますが、
真ん中あたりで通常とは違った動きをしているのがわかります。
今回は特異スペクトル変換法を使用して、この異常部分を発見したいと思います。

特異スペクトル変換法


特異スペクトル変換法は以下の図のように説明されます。

(井手剛・杉山将(2015)「異常検知と変化検知(機械学習プロフェッショナルシリーズ)」,講談社.を参考に著者作成)

特異スペクトル変換法は、ある時刻tの周りに部分時系列を現在側と過去側に作成します。
現在側の部分時系列をテスト行列(Z^(t))、過去側の部分時系列を履歴行列(X^(t))と呼びます。
このテスト行列、履歴行列にはそれぞれある一定区間に区切られたWindowが複数含まれています。

時刻tの変化度(異常度)はテスト行列に含まれるWindow群の波形と
履歴行列に含まれるWindow群の波形を比較して、どれだけ食い違っているかで定義されます。
ただし、履歴行列とテスト行列を比較する前にそれぞれのWindow群の特徴を表す代表的な波形を数本作成しておきます。

この代表的な波形を作成する作業が特異値分解(SVD)になります。
したがって、履歴行列、テスト行列それぞれを代表する特徴的な波形数本同士を比較し、
どれほど食い違いがあるかを調べます。

特異値分解(Singular Value Decomposition)


特異値分解は固有値分解を正方行列だけでなく、長方形行列にも適用できるように拡張したものと考えられます。
今回は履歴行列、テスト行列に含まれるWindowを列に並べて、
全体の情報をある程度保ちながら次元削減するために使用します。

例えば履歴行列、テスト行列にそれぞれ50本の波形が含まれていたとして(すなわち50Window)、
波形1本が1列に対応するように並べた後に、特異値分解を使用して次元削減します。
イメージとしてはPCAと同じような働きをすると考えるとわかりやすいです。
50本の波形全体の情報を保ちつつ、次元を減らすわけです。
この操作を先ほど「Window群の特徴を表す代表的な波形を数本作成する」と表現しました。

特異値分解は元データのXを上記のように3つの部分に分解します。
式で書くと以下ようになります。

$$ X = U\Gamma^{1/2}V^T$$

Uの各列を左特異ベクトル、Gamma^{1/2}の対角成分を特異値、Vの各列を右特異ベクトルと呼びます。
左特異ベクトルの行列から、いくつか列を抜き出す事によって、次元削減を行います。
履歴行列、テスト行列を特異値分解し、左特異ベクトルからそれぞれr列、m列抜き出したものを
$$U_r^{(t)}, Q_m^{(t)}$$とします。

異常度は以下のように定義されます。
$$a(t) = 1 – (U_r^{(t)} Q_m^{(t)}の最大特異値)$$

RapidMinerでの実装


特異スペクトル変換法をRapidMinerとpythonで実装した例が以下になります。
ChangeDetection_SST.rmpというファイルをRapidMinerで開いてください。
(メニューバー>File>Import Processより上記rmpファイルを選択してください。)

以下にプロセスを解説していきます。

①Retrieve ECG_include_anomaly_3000
データを読み込むオペレータです。
プロセスを保存した場所と同じ場所にECG_include_anomaly_3000を取り込んでください。

②Macro NumOfExample(Extract Macro)
データセットをのサンプル数をカウントして格納しておくマクロです。
ここではNumOfExampleに3000が格納されます。
マクロについては以下のブログを参考にしてください。
あなたの前処理を早くする!前処理の三種の神器【マクロ編】

③Windowing
Window sizeを100, step sizeを1として、Windowを作成します。
1Windowにつき1行、Window内のデータをそのまま列に並べます。
サンプル数はN – Window size + 1すなわち3000 – 100 + 1 = 2901
になります。
Windowingについては以下のブログも参考にしてください。
Time Series as Supervised Learning: Lag Variables

④Make Macros(Subprocess)
中にExtract MacroとGenerate Macroが配置されています。
いずれもマクロを作成するためのものです。

ここでは以下のマクロを作成しています。
・window_size : 窓の大きさ。ここでは100が格納されます。
・k : 履歴行列とテスト行列に含まれるWindowの数。ここでは50が格納されます。
・L : ラグ。履歴行列とテスト行列をどれだけずらすか。ここでは25が格納されます。
・iternum : (3000 – 100 + 1) – (50 + 25) + 1 = 2827回
・m : 特異値分解後、左特異ベクトルの行列から何列使用するか。3を格納しています。

⑤Loop
ここで特異値分解を計算しています。
中は以下のようになっています。

Generate Macroで履歴行列のスタート地点と終了地点、テスト行列のスタート地点と終了地点を設定しています。反復ごとに1時点ずれていきます。
Transposeオペレータを置いているのは、1列1windowにするためです。
Transpose後、履歴行列もテスト行列も100行50列のデータになります。
SVDオペレータはデータを特異値分解し、左特異ベクトルの行列を指定した列数で返却します。
列数はdimensionality reductionパラメータとdimensionパラメータで設定しています。
Generate ID(Generate Attribute)オペレータでそれぞれのSVDの結果にID番号(反復回数)を割り当てます。

⑥Append
Loopの結果を行結合して、一つにまとめます。

⑦Execute Python
ID番号が一致する履歴行列とテスト行列の内積を計算し、その最大特異値を1から引いたものを異常度と定義します。
コードは以下の通りです。

⑧Store
異常度をanomaly_scoreという名前で保存します。

元の波形と異常度の関係はcompare_anomaly_and_seriese.rmpで確認できます。
異常な動きをしている所の異常度が著しく大きくなっています。

いかがでしたでしょうか。
皆さんも一度、RapidMinerを使用して異常検知や変化検知を試してみてください。

=========================================

▼1日で機械学習の基礎が学べる!無償ハンズオンセミナー随時開催中▼
URL: https://www.ksk-anl.com/event/
▼RapidMinerをもっと知りたい方はこちら▼
URL: https://www.rapidminer.jp/
▼早速使ってみたい!RapidMinerの体験版の無償DL▼
URL: https://www.rapidminer.jp/download/
▼導入を検討中!わかりやすい詳細なQ&A▼
URL: https://www.rapidminer.jp/faq/

==========================================

$$ \begin{aligned} \newcommand\argmin{\mathop{\rm arg~min}\limits} \boldsymbol{\beta}_{\text{ridge}} & = \argmin_{\boldsymbol{\beta} \in \mathcal{R^p}} \biggl[ ||\boldsymbol{y}-\boldsymbol{X\beta}||^2 + \lambda ||\boldsymbol{\beta}||^2 \biggr] \\ & = (\boldsymbol{X}^T\boldsymbol{X} + \lambda\boldsymbol{I_{p+1}})^{-1}\boldsymbol{X}^T\boldsymbol{y} \end{aligned} $$
PAGE TOP