概要
先のGaussianMixture modelに続き、異常検知(Anomaly detection)について調べていて発見した副産物について残す。
今回試したのはマハラノビス距離をmetricsとするMinimum Covariance Determinant estimator (以下MCD)による異常検知。
クラスタリングの話ではない。
- 実施期間: 2021年8月
- 環境：Python3.8
- パケージ：scikit-learn
MCDとマハラノビス距離
Anomaly detectionについて調査中、scikit learnのMCDチュートリアルを見つけた。
以前、scipyのマハラノビスを試したときはいまひとつ２つの集まりを分離できなかったので、ダメだこりゃ、と思っていたがMinimum Covariance Determinant estimator(MCD)なる手法で格段に分離性能が向上することを知った。
汚染されたoutlierが全dataの半数近く($\frac{n_{samples}-n_{features} - 1}{2}$)あってもその反対の半分近い($\frac{n_{samples}-n_{features} + 1}{2}$)のdatasetで最小の共分散を持つdatasetを探せれば、outlierを区別できるよね、というロバストなモデルらしい。
感覚的には、ひとつの集合に見えるふたつの異なる集合間で、どれほど相関があるか共分散を使って分離させる。そうすると両集合内ではdate間に相関があり、異なる集合間では相関がなくなる。date単体で評価したいのでマハラノビス距離が登場する。つまり相関の強い集合内のdata点間はその距離が短くなり相関が弱い集合外のdata点の距離は遠くなる。距離が遠い＝特異点、ってな感じ。
準備
scikit learnの例題は目的が異なりやや冗長なので、要点を絞ってやってみる。
datasetはみんな大好きアヤメちゃんから一部拝借した。下図の緑全部(50点)と一部かぶった青の無作為で一部(10点)を使用する。
df_MCD = df_iris[df_iris.target != 0.0] # Target1.0, 2.0
df_MCD = df_MCD.drop(range(110,150)) # Target 2.0の40点はomit
# 元dataは4次元だが、図示するために2次元とする。
X = np.array(df_MCD[['petal width (cm)','petal length (cm)']])
# x_train = np.array(df_iris[['petal width (cm)','petal length (cm)', 'sepal length (cm)', 'sepal width (cm)']])
Y = np.array(df_MCD['target'])
color_codes = {0:'r', 1:'g', 2:'b'}
colors_pred = [color_codes[x] for x in Y]
plt.xlabel('width', fontsize=12)
plt.ylabel('length', fontsize=12)
plt.scatter(X[:,0], X[:,1], c=colors_pred)
一部が重なり、かつ緑・青の両方でwidth - lengthに同じ相関があるという、オイラがアヤメちゃんで一番お気に入りの部分。色がついてなければ同じ集合だろうと思ってしまう。
MCDによるAnomaly Detection
早速準備した全60点のdatasetをMCDに入力し、その共分散行列を確認する。
import matplotlib.pyplot as plt
from sklearn.covariance import MinCovDet
# fit a MCD robust estimator to data
robust_cov = MinCovDet().fit(X)
print('Estimated covariance matrix:\n', robust_cov.covariance_)
Estimated covariance matrix:
[[0.04327566 0.07769319]
[0.07769319 0.21426374]]
共分散行列はマジョリティである緑の集合のものとなっているはずである。そこで緑と青の共分散行列をもとめる。
print(df_MCD.loc[50:100,['petal width (cm)','petal length (cm)']].cov())
print(df_MCD.loc[10:110,['petal width (cm)','petal length (cm)']].cov())
つまり、おおよそ緑を特定しているといえよう。
マハラノビス距離でMCDの結果をみてみる。
fig, ax = plt.subplots(figsize=(10, 5))
# Plot data set
inlier_plot = ax.scatter(X[:, 0], X[:, 1],
color='blue', label='inliers')
outlier_plot = ax.scatter(X[:, 0][-10:], X[:, 1][-10:],
color='red', label='outliers')
ax.set_xlim(ax.get_xlim()[0], ax.get_ylim()[0])
ax.set_title("Mahalanobis distances of a contaminated data set")
# Create meshgrid of feature 1 and feature 2 values
xx, yy = np.meshgrid(np.linspace(plt.xlim()[0], plt.xlim()[1], 100),
np.linspace(plt.ylim()[0], plt.ylim()[1], 100))
zz = np.c_[xx.ravel(), yy.ravel()]
# Calculate the MCD based Mahalanobis distances
mahal_robust_cov = robust_cov.mahalanobis(zz)
mahal_robust_cov = mahal_robust_cov.reshape(xx.shape)
robust_contour = ax.contour(xx, yy, np.sqrt(mahal_robust_cov),
levels=np.logspace(0, 3, 20),
linestyles='dashed', linewidths=0.5)
robust_contour.clabel(fmt='%1.1f', fontsize=12)
plt.xlabel('width', fontsize=12)
plt.ylabel('lenght', fontsize=12)
# Add legend
ax.legend([robust_contour.collections[1],
inlier_plot, outlier_plot],
['MCD dist', 'inliers', 'outliers'],
loc="upper right", borderaxespad=0)
plt.show()
マハラノビス距離で見ると良い感じに区別できているが、実際はやはりdata見ながら閾値次第ってところか。
重なったところは無理があるが、それは50点+10点ということを知ってて色をつけているからであって、実際にはそういった情報はないから、そうするとdataの背景にある平均・分散を推定し識別してくれるということには大きな有用性がある。
もちろん各data点の距離も計算してくれる。
mahal_robust_cov = robust_cov.mahalanobis(X)
print(mahal_robust_cov)
[ 0.94672455 0.62971656 1.40848119 0.76069509 0.45462248 0.91842642
1.44516541 4.34352419 1.68954317 3.49310377 3.05434983 2.63413978
4.14557506 0.94672455 5.07193291 0.05269056 0.62971656 5.10339059
0.62971656 1.43930893 6.20138769 0.76069509 1.40848119 5.80742846
0.1157634 0.05269056 1.73778287 2.55886182 0.62971656 3.05434983
1.4476458 2.75126944 0.78412102 2.30369632 0.62971656 2.49504092
0.52605189 0.39383316 0.29919436 0.76069509 2.07476196 0.40218973
0.54920222 4.34352419 0.08421714 0.8189351 0.08421714 0.1157634
10.38918659 0.29919436 31.86608188 7.42248905 12.15992136 6.2653267
15.72251362 20.40013106 5.70015147 19.26375677 8.7465464 30.65778604]
以上