はじめに
本連載について
こんにちは，(株)日立製作所 Lumada Data Science Lab. の露木です。
化学プラントや工場設備，あるいはもっと身近なモーターや冷蔵庫などの故障予兆検知を行う際に，振動や音，温度，圧力，電圧，消費電力のような値を取得できる複数のセンサーで測定した多次元の時系列数値データを分析することがあります。このような故障予兆検知は，機械学習の分野では異常検知問題として解くことができます。
そこで本連載では時系列数値データの異常検知を最終目標とし，数回の記事に分けてアルゴリズムの基礎的な説明と実装を示していきます。初回の本稿では，後に説明する異常検知アルゴリズムの出発点となるマハラノビス距離について取り上げます。
- 第1回. scipyを使って特徴量の相関を考慮したマハラノビス距離を計算する (本稿)
- 第2回. ホテリングの$T^2$法による多変量正規分布を仮定した異常検知
- 第3回. GMMによる多峰性分布にもとづいたデータの異常検知
- 第4回. 移動窓を使った多次元時系列データの異常検知
想定シナリオ
まず最初に何かの機器が故障する具体的なシナリオを想定して，時系列数値データの異常検知問題に必要な要素技術を考えます。
本連載では，図1のように2つの熱源と冷却ファン，2つの温度センサがあるような仮想的な機器を想定します。この機器では，熱源1に接続された運転スイッチが時刻とともに運転，停止，運転，停止と繰り返しており，それに対応して熱源1の温度も上下していきます。一方，熱源2は常に運転状態にありますが，冷却ファンによって一定以下の温度に抑えられています。
ここで，熱源2の冷却ファンが故障したとします。この故障を2つの温度センサの値から見つけ出すことが最終目標です。
図1. 想定シナリオ
2つの温度センサの値を，横軸に時刻をとった時系列数値データとして可視化すると図2のようなグラフになります。正常時は熱源1の運転状態に応じて，2つの温度センサの値が連動して上下します。そして時刻150で冷却ファンが故障すると，熱源2の加熱に対応して温度センサ2の値のみが上昇する異常値が記録されます。
図2. 時系列データの例
時系列データを各時刻でのセンサ値の組を持つベクトルとして捉え，散布図の形式で可視化すると，異常値の存在がより明確になります。図2と同じデータについて，今度は縦軸・横軸に2つの温度センサの値をとって可視化すると図3のような散布図を得られます。このように可視化すると，データ分析によって熱源1の「運転モード」と「停止モード」，そしてその間の「移行期」および熱源2が過熱した「故障モード」の4つに分類できるように見えてきます。
図3. 期待される分析結果
利用可能なデータは温度センサの値のみであるため，機械学習の言葉で言えば教師なしの分類問題を解くことになります。より具体的に言えば，クラスタリングによって与えられた特徴量空間における「正常なデータ」の分布を学習し，そこから遠く離れたデータを「異常なデータ」とみなせば良いと言えます。これを実現するためには要素技術として以下の4つが必要です。
- 特徴量空間における距離の定義
- 正常データの分布を仮定し，異常なデータと正常なデータを分ける閾値を設定する方法
- 複数のモードに対応できるクラスタリング
- 時系列データに対応するための移動窓の活用
上記の1～4を実装する方法について，先に示した本連載の第1回～第4回で順に説明していきます。
今回扱う内容
初回である本稿では特徴量空間における距離の定義を扱います。距離の定義は，正常なデータと異常なデータを正しく分離するために重要な役割を果たします。
図3のように横に細長い「移行期」や「故障モード」を分離したクラスタリングは，等方的なユークリッド距離では難しく，マハラノビス距離が有効です。マハラノビス距離は母集団の分布を考慮し，特徴量の間に強い相関がある場合でも正常データと異常データを分離しやすい特徴を持った距離の定義方法です。実際に，マハラノビス距離は異常検知のアルゴリズムとして有名なホテリングの $T^2$ 法やGMM (混合正規分布モデル) の内部で使われています。その理論的な詳細は成書 1 や多くのWeb記事 2 3 4 に譲るとして，本稿ではscipyの距離関数 5 を利用した実装方法を示します。
マハラノビス距離の特徴
ユークリッド距離の欠点
距離の定義としてユークリッド距離はもっとも著名であり，比較対象として適任といえます。${\bf x}$と${\bf y}$の2つのベクトル間のユークリッド距離$d_e$は次式で計算できます。
d_{e} = \sqrt{ \sum_i ({x}_i - {y}_i)^2}
しかし，先に述べたようにユークリッド距離は背景にある母集団の分布を考慮した距離ではありません。特に，特徴量に強い相関がある場合，ユークリッド距離にもとづいたクラスタリングではうまく分離できないことがあります。例えば図4のように，原点 (青色のデータ群の標本平均) を中心としたユークリッド距離による分離を考えます。図4では，ユークリッド距離 $d_e=5$ の等高線では青色のデータ群を包み込むには小さすぎる一方で， $d_e=10$ ではオレンジ色のデータ群も含んでしまいます。
図4. ユークリッド距離で分離する場合
マハラノビス距離の利点
このような問題に対処するため，異常検知の分野ではマハラノビス距離がよく使われます。 1 マハラノビス距離$d_m$は分散共分散行列$\Sigma$を用いて次式で定義されます。
d_{m} = \sqrt{ ({\bf x} - {\bf y})^{\rm T} \Sigma^{-1} ({\bf x} - {\bf y}) }
定義式に母集団の分散が含まれることから感覚的にも理解できますが，以下の図5のように原点を中心としたマハラノビス距離の等高線図は分散の大きい方向に引き伸ばされた形式になります。これにより，ユークリッド距離では難しかった青色のデータ群とオレンジ色のデータ群の分離が可能になります。
図5. マハラノビス距離で分離する場合
scipy (Python) による実装
マハラノビス距離の計算方法
最初に結論を述べると，scipyに組み込みの関数 scipy.spatial.distance.mahalanobis() を使えば，以下のように簡単にマハラノビス距離を計算できます。
ただし，numpyのcov関数 はデフォルトで不偏分散を計算する (つまり，$1/(N-1)$で行列要素が規格化されている) ことに注意してください。cov関数に bias=True の引数を与えれば標本分散を計算します (つまり，$1/N$で規格化する)。 6
import numpy as np
from scipy.spatial import distance
X = np.array([[0, 0], [1, 1], [2, 2]])
# 分散共分散行列を計算
cov = np.cov(X.T)
# 分散共分散行列の逆行列を計算
cov_i = np.linalg.pinv(cov)
# 2つの標本 [1, 1] と [0, 0] のマハラノビス距離を計算する
d = distance.mahalanobis([1, 1], [0,0], cov_i)
print("マハラノビス距離の計算結果: %1.2f" % d)
マハラノビス距離の計算結果: 1.00
なお，scipyにはユークリッド距離を計算する関数 scipy.spatial.distance.euclidean() も用意されており，こちらも簡単に計算できます。 2つの標本 [0, 0] と [1, 1] のマハラノビス距離は 1.00 でしたが， ユークリッド距離は 1.41 となり，母集団の分散に応じてマハラノビス距離で数えたほうがユークリッド距離よりも近いということになりました。
# 2つの標本 [1, 1] と [0, 0] のユークリッド距離を計算する
d = distance.euclidean([1, 1], [0, 0])
print("ユークリッド距離の計算結果: %1.2f" % d)
ユークリッド距離の計算結果: 1.41
トイデータを使った2種類の距離の比較
トイデータの性質
マハラノビス距離の計算方法を示したので，次に本稿の冒頭で示した等高線図を描いてユークリッド距離とマハラノビス距離の比較を行います。まずは人工的な右肩上がりの成分に正規分布のノイズを足したトイデータを作成します。
X1 = np.array([(x, x + np.random.normal()) for x in np.arange(-5, 5, 0.1)])
X2 = np.array([(x, 10 + x + np.random.normal()) for x in np.arange(-10, 0, 0.1)])
このトイデータは以下の図6のように右肩上がりに偏った分布のデータセットです。
import matplotlib.pyplot as plt
# グラフ設定
myfig = plt.figure()
myax = myfig.add_subplot(1, 1, 1)
myax.set_aspect('equal')
plt.ylim(-20, 20)
plt.xlim(-20, 20)
# 散布図を描画
plt.plot(X1[:,0], X1[:,1], '.')
plt.plot(X2[:,0], X2[:,1], '.')
plt.show()
図6. データセットの分布
ユークリッド距離では適切に分離できない
まず，青色のデータ群の中心点 (標本平均) からのユークリッド距離を等高線図として可視化します。図7を見ると，やはりユークリッド距離の等高線は2つのデータ群にまたがっており，回のトイデータのような強い相関を持った特徴量空間に広がった2つのデータ群をうまく分離できないといえます。
# 青色のデータ群の中心 (標本平均) を計算
mean = X1[:,0].mean(), X1[:,1].mean()
# z軸値 (標本平均からのユークリッド距離) を計算
z = np.array(
[
[(i, j, distance.euclidean([i, j], mean))
for i in np.linspace(-20, 20, 100)]
for j in np.linspace(-20, 20, 100)]
)
# グラフ設定
myfig = plt.figure()
myax = myfig.add_subplot(1, 1, 1)
myax.set_aspect('equal')
plt.ylim(-20, 20)
plt.xlim(-20, 20)
# 等高線図を描画
cont = plt.contour(z.transpose()[0], z.transpose()[1], z.transpose()[2], levels=[5,10,15], colors=['k'])
cont.clabel(fmt=' $d_e$=%1.1f ', fontsize=12)
# 散布図を描画
plt.plot(X1[:,0], X1[:,1], '.')
plt.plot(X2[:,0], X2[:,1], '.')
plt.show()
図7. ユークリッド距離の計算結果
マハラノビス距離なら適切に分離できる
次に同じ可視化をマハラノビス距離について行います。可視化で得た図8を見れば，マハラノビス距離であれば母集団の分布を考慮するので$d_m=5$の等高線で2つのデータ群を適切に分離できていることがわかります。
# 分散共分散行列を計算
cov = np.cov(X1.T)
# z軸値 (標本平均からのマハラノビス距離) を計算
z = np.array(
[
[(i, j, distance.mahalanobis([i,j], mean, np.linalg.pinv(cov)))
for i in np.linspace(-20, 20, 100)]
for j in np.linspace(-20, 20, 100)]
)
# グラフ描画
myfig = plt.figure()
myax = myfig.add_subplot(1, 1, 1)
myax.set_aspect('equal')
plt.ylim(-20, 20)
plt.xlim(-20, 20)
# 等高線図を描画
cont = plt.contour(z.transpose()[0], z.transpose()[1], z.transpose()[2], levels=[5,10,15], colors=['k'])
cont.clabel(fmt='$d_m$=%1.1f', fontsize=12)
# 散布図の描画
plt.plot(X1[:,0], X1[:,1], '.')
plt.plot(X2[:,0], X2[:,1], '.')
plt.show()
図8. マハラノビス距離の計算結果
おわりに
本稿で計算したマハラノビス距離を使えばデータを分離できますが，これだけではデータの中心や異常データと正常データを分ける閾値を決定できません。そこで，次回の記事では，ホテリングの$T^2$法 1 を実装する方法を説明します。ホテリングの$T^2$法は，正常データについて正規分布を仮定した上で，データが出現する確率にもとづいて異常度の閾値を決定して異常検知を行う方法です。なお，ホテリングの$T^2$法では異常度としてマハラノビス距離の2乗値を使います。
また，混合正規分布モデル (GMM) によるクラスタリングでも一般にマハラノビス距離が使われます。一方，クラスタリングの手法として著名なK-Means法は，ユークリッド距離を内部的に使う事が多いアルゴリズムです。一般に，K-meansよりもGMMのほうが性能が良くなることが多い 7 のですが，その一因としてマハラノビス距離の採用が挙げられます。