#はじめに
この記事ではマハラノビス距離の計算と可視化を通して、
マハラノビス距離を用いた異常検知のイメージを確認します。
理論や数式は省いて、ソースコードと可視化を中心に、定性的な理解に重きを置いています。

この記事ですること
・相関のあるデータの生成と可視化
・分散共分散行列の計算と可視化
・マハラノビス距離の計算と可視化

マハラノビス距離の定義等については、既にQiita内でも解説してくださっている方がいらっしゃるので、
そちらを参考にされると良いと思います。
[教師なし学習による異常値検知: マハラノビス距離 (理論編)]
(https://qiita.com/shopetan/items/ceb7744facc21c3881d2)
[scipyを使って特徴量の相関を考慮したマハラノビス距離を計算する]
(https://qiita.com/MasafumiTsuyuki/items/de19d8ec274e961ec946)

scikit-learnでの公式ページの解説はこちらからご確認ください。
[Robust covariance estimation and Mahalanobis distances relevance]
(https://scikit-learn.org/stable/auto_examples/covariance/plot_mahalanobis_distances.html)

本論

ライブラリ

バージョン確認

import numpy
import matplotlib
import scipy
import sklearn
print("numpy: {0}".format(numpy.__version__))
print("matplotlib: {0}".format(matplotlib.__version__))
print("scipy: {0}".format(scipy.__version__))
print("sklearn: {0}".format(sklearn.__version__))

筆者環境では以下のバージョンがインストールされております。

numpy: 1.17.0
matplotlib: 3.2.2
scipy: 1.5.0
sklearn: 0.23.1

使用する部分のインポート

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from sklearn.covariance import EmpiricalCovariance

相関のあるデータの生成と可視化

今回は正の相関を持つ2次元正規分布のデータを生成して学習データとして利用します。
分散をそれぞれ1.0、共分散を0.7に設定しています。

mean = np.array([0,0])  # 平均
cov_org = np.array([[1.0,0.7],[0.7,1.0]])  # 分散共分散行列  
np.random.seed(0)
X = multivariate_normal(mean, cov_org).rvs(size=1000)  # データ生成

生成したデータを散布図で可視化します。

fig,ax = plt.subplots(figsize=(8,6))
ax.plot(X[:,0], X[:,1] ,'.')
ax.set_xlim([-4,4])
ax.set_ylim([-4,4])
ax.set_xlabel('X[0]')
ax.set_ylabel('X[1]')
ax.set_aspect('equal')
ax.set_title('Fig.1')
plt.savefig('fig1.png')

正の相関を持つデータが生成できていることが確認できます。

分散共分散行列の計算と可視化

生成したデータから分散共分散行列を再計算してみます。
これがデータの「学習」に相当します。

empCov = EmpiricalCovariance().fit(X)
cov_cal = empCov.covariance_

ヒートマップで可視化してみます。

fig, ax = plt.subplots(figsize=(8,6))
im = ax.imshow(cov_cal, vmin=0, cmap='Reds')
cbar = ax.figure.colorbar(im, ax=ax)
ax.set_xticks(np.arange(cov_cal.shape[0]))
ax.set_yticks(np.arange(cov_cal.shape[1]))
ax.tick_params(top=True, bottom=False,labeltop=True, labelbottom=False)

for i in range(cov_cal.shape[0]):
    for j in range(cov_cal.shape[1]):
        text = ax.text(i, j, cov_cal[i, j], ha='center', va='center', color='w')
ax.set_title('Fig.2')
plt.savefig('fig2.png')

データの生成時に設定した分散共分散行列に従ってデータが生成されていることが確認できます。

マハラノビス距離の計算と可視化

テスト用にデータを生成します。
今回は[-2,-2]～[2,2]の範囲で1間隔の格子状にデータを生成しています。
説明の都合上、A:[-2,-2],B:[-1,-2],C:[0,-2],…,M:[0,0],…,Y:[2,2]と対応付けることにします。

Y0,Y1 = np.meshgrid(np.arange(-2,3),np.arange(-2,3))
Y = np.array([Y0.reshape(-1),Y1.reshape(-1)]).T
labels = ['A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y']

これを先ほど使用したEmpiricalCovarianceクラスのmahalanobis()に入れるとマハラノビス距離の2乗が計算できます。

md = np.sqrt(empCov.mahalanobis(Y))
print('Mahalanobis Distance')
print(md)

実行結果

Mahalanobis Distance
[2.24902063 2.12256063 2.84114155 3.96970512 5.25020799 2.16386194
 1.13576907 1.42031549 2.61997568 3.97880646 2.88870249 1.45353983
 0.0246127  1.41706927 2.85223012 4.00562162 2.64077329 1.42197552
 1.09086983 2.11604144 5.27100616 3.98405422 2.84280177 2.10000179
 2.20412   ]

棒グラフで可視化するとこのようになります。

fig,ax = plt.subplots(figsize=(8,6))
ax.bar(labels, md)
ax.set_xlabel('Label')
ax.set_ylabel('Mahalanobis Distance')
ax.set_title('Fig.3')
plt.savefig('fig3.png')

この結果を詳しく見ていくことにします。
先ほどの散布図に、テスト用データおよびラベルと、マハラノビス距離の等高線も含めて可視化します。

### 等高線用データ
cX, cY = np.meshgrid(np.linspace(-4, 4, 100), np.linspace(-4.0, 4.0, 100))
cZ = np.sqrt(empCov.mahalanobis(np.array([cX.reshape(-1),cY.reshape(-1)]).T)).reshape(100,100)
### 散布図 + 等高線
fig,ax = plt.subplots(figsize=(8,6))
ax.plot(X[:,0], X[:,1], '.', alpha=0.3)  # 学習データ
ax.plot(Y[:,0], Y[:,1], 'x', color='red')  # テストデータ
ax.clabel(ax.contour(cX, cY, cZ, np.arange(0, 10, 1), cmap='jet'), inline=True, fontsize=10)  # 等高線用データ
for i in range(Y.shape[0]):
    text = ax.text(Y[i,0]+0.1, Y[i,1]+0.1, labels[i], color='r')
ax.set_xlim([-4,4])
ax.set_ylim([-4,4])
ax.set_xlabel('X[0]')
ax.set_ylabel('X[1]')
ax.set_aspect('equal')
ax.set_title('Fig.4')
plt.savefig('fig4.png')

M[0,0]の点は分布の中央に位置しているので、マハラノビス距離はほぼ0となります。
つまり、異常度は極めて低いと言うことができます。
E[2,-2]とU[2,-2]が相関から外れているため、マハラノビス距離が5.3程度と大きくなっています。
そのため、例えば「マハラノビス距離5以上を異常とみなす」と設定した場合は、
E・Uの2点を異常として検知することができます。
また、A[-2,-2]とE[2,-2]を比べてみると、中心からの画像上の距離(=ユークリッド距離)は同じですが、
Aの方は相関から外れていないのでマハラノビス距離も2.2と小さくなっています。
他に、A[-2,-2]とI[1,-1]を比べてみると、ユークリッド距離が小さいのはIですが、相関から外れているため、
マハラノビス距離は2.6とAより大きくなっています。
以上のように相関を考慮した異常検知ができていることが確認できると思います。

まとめ

以上、コードと可視化を交えたマハラノビス距離の解説でした。
どなたかのお役に立てれば幸いです。
お気づきの点がありましたらご指摘ください。

付録

使用したコード全体

## ライブラリ
### バージョン確認
import numpy
import matplotlib
import scipy
import sklearn
print("numpy: {0}".format(numpy.__version__))
print("matplotlib: {0}".format(matplotlib.__version__))
print("scipy: {0}".format(scipy.__version__))
print("sklearn: {0}".format(sklearn.__version__))

### 使用する部分のインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from sklearn.covariance import EmpiricalCovariance


## データの生成と可視化
### データの生成
mean = np.array([0,0])  # 平均
cov_org = np.array([[1.0,0.7],[0.7,1.0]])  # 分散共分散行列  
np.random.seed(0)
X = multivariate_normal(mean, cov_org).rvs(size=1000)  # データ生成

###　散布図
fig,ax = plt.subplots(figsize=(8,6))
ax.plot(X[:,0], X[:,1] ,'.')
ax.set_xlim([-4,4])
ax.set_ylim([-4,4])
ax.set_xlabel('X[0]')
ax.set_ylabel('X[1]')
ax.set_aspect('equal')
ax.set_title('Fig.1')
plt.savefig('fig1.png')


## 分散共分散行列の計算と可視化
### 分散共分散行列の再計算(学習)
empCov = EmpiricalCovariance().fit(X)
cov_cal = empCov.covariance_

### ヒートマップ
fig, ax = plt.subplots(figsize=(8,6))
im = ax.imshow(cov_cal, vmin=0, cmap='Reds')
cbar = ax.figure.colorbar(im, ax=ax)
ax.set_xticks(np.arange(cov_cal.shape[0]))
ax.set_yticks(np.arange(cov_cal.shape[1]))
ax.tick_params(top=True, bottom=False,labeltop=True, labelbottom=False)

for i in range(cov_cal.shape[0]):
    for j in range(cov_cal.shape[1]):
        text = ax.text(i, j, cov_cal[i, j], ha='center', va='center', color='w')
ax.set_title('Fig.2')
plt.savefig('fig2.png')


## マハラノビス距離の計算と可視化
### テスト用データ生成(格子点)
Y0,Y1 = np.meshgrid(np.arange(-2,3),np.arange(-2,3))
Y = np.array([Y0.reshape(-1),Y1.reshape(-1)]).T
labels = ['A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y']

### マハラノビス距離計算
md = np.sqrt(empCov.mahalanobis(Y))
print('Mahalanobis Distance')
print(md)

### 棒グラフ
fig,ax = plt.subplots(figsize=(8,6))
ax.bar(labels, md)
ax.set_xlabel('Label')
ax.set_ylabel('Mahalanobis Distance')
ax.set_title('Fig.3')
plt.savefig('fig3.png')

### 等高線用データ
cX, cY = np.meshgrid(np.linspace(-4, 4, 100), np.linspace(-4.0, 4.0, 100))
cZ = np.sqrt(empCov.mahalanobis(np.array([cX.reshape(-1),cY.reshape(-1)]).T)).reshape(100,100)
### 散布図 + 等高線
fig,ax = plt.subplots(figsize=(8,6))
ax.plot(X[:,0], X[:,1], '.', alpha=0.3)  # 学習データ
ax.plot(Y[:,0], Y[:,1], 'x', color='red')  # テストデータ
ax.clabel(ax.contour(cX, cY, cZ, np.arange(0, 10, 1), cmap='jet'), inline=True, fontsize=10)  # 等高線用データ
for i in range(Y.shape[0]):
    text = ax.text(Y[i,0]+0.1, Y[i,1]+0.1, labels[i], color='r')
ax.set_xlim([-4,4])
ax.set_ylim([-4,4])
ax.set_xlabel('X[0]')
ax.set_ylabel('X[1]')
ax.set_aspect('equal')
ax.set_title('Fig.4')
plt.savefig('fig4.png')