23
23

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

GaussianMixture modelによる異常検知 備忘録

Last updated at Posted at 2021-06-02

#概要
異常検知(Anomaly detection)について調べていて発見した副産物について書き残そう。
結局オーソドックスなVAEでいくことにしたのだが、このGMMの方がセンスが良く感じる。
DNNを使ったときの手抜き感は小生だけだろうか?

  • 実施期間: 2021年5月
  • 環境:Google Colaboratory
  • パケージ:scikit-learn

なお、日本語で「異常」というとabnormalのイメージがあるが、abnormalとanomalyはニュアンスが違うので、職場では「異常」という言葉は使わないようにしている。

【変更履歴】
2021/9/4 尤度スコアで等高線を追加した。

##モチベーション
Anomaly detectionについて調査中、GaussianMixture Model(GMM)についてのブログを見つけた。
クラスタリングの一種だがK-MeansのようにDeterministicに距離でクラスタ分けするようなアルゴリズムではなく、各点がどのようなGaussian分布に属すべきかをEM(Expectation maximization)で最尤推定するものらしい。
教科書にしているO'REILLYのHands-onでも詳しく説明されており、面白かったので少し味見してみることにした。

以下の順に備忘録として残す。

  1. 準備
  2. K-Meansによるクラスタリングの復習
  3. GMMによるクラスタリング
  4. GMMによるクラスタリング数のもとめ方
  5. GMMによるAnomaly Detection

##1. 準備
いつもimportしているパケージから取捨選択

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn import metrics
from sklearn import linear_model
from sklearn.metrics import plot_confusion_matrix
from IPython.core.debugger import Pdb # 念のため

今回使用するdatasetはみんな大好きアヤメちゃん

iris = load_iris()
feature_names = iris.feature_names
target_names = iris.target_names

df_iris = pd.DataFrame(data = np.c_[iris['data'], iris['target']],
                     columns = iris['feature_names'] + ['target'])

X = df_iris[['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']]
y = df_iris['target']

x_train = np.array(df_iris[['petal width (cm)','petal length (cm)']])
#x_train = np.array(df_iris[['sepal length (cm)', 'sepal width (cm)']])
y_test = np.array(df_iris['target'])

color_codes = {0:'r', 1:'g', 2:'b'}
colors_pred = [color_codes[x] for x in y_test]
plt.scatter(x_train[:,0], x_train[:,1], c=colors_pred)

image.png

また、評価結果用の関数を定義

def print_evaluate(x_train, y_test, y_pred, ave='macro'):
    print('accuracy = \t', metrics.accuracy_score(y_true=y_test, y_pred=y_pred))
    print('precision = \t', metrics.precision_score(y_true=y_test, y_pred=y_pred, average=ave))
    print('recall = \t', metrics.recall_score(y_true=y_test, y_pred=y_pred, average=ave))
    print('f1 score = \t', metrics.f1_score(y_true=y_test, y_pred=y_pred, average=ave))
    print('confusion matrix = \n', metrics.confusion_matrix(y_true=y_test, y_pred=y_pred))

※2値のときは average='binary'

##2. K-Meansによるクラスタリング
色々とアルゴリズムはあるけれど、GMMで少し言及するK-Meansについて復習する。
簡単に言えば点群の集まりの中心(centroid)がどこかを探すやつ。
指定数のcentroidと全sampleの距離が小さくなるように全centroidを移動させ、収束したときの各centroidに近いsampleがそのクラスタに属するとするアルゴリズム
まず、オーソドックスな**K-Means**を見てみる。クラスタ(centroid)数は既知のn_clusters=3とする。

from sklearn.cluster import KMeans
from sklearn.cluster import MiniBatchKMeans

## K-Means
kmeans = KMeans(n_clusters=3, random_state=0).fit(x_train)
y_pred = kmeans.labels_
y_cent = kmeans.cluster_centers_

# 精度評価結果
print_evaluate(x_train, y_test, y_pred, 'macro')

colors_pred = [color_codes[x] for x in y_pred]
plt.scatter(x_train[:,0], x_train[:,1], c=colors_pred)
plt.scatter(y_cent[:,0], y_cent[:,1], s = 250, marker='*',c='y')
plt.show

accuracy = 0.37333333333333335
precision = 0.37286324786324787
recall = 0.37333333333333335
f1 score = 0.3730825663598773
confusion matrix =
[[50 0 0]
[ 0 2 48]
[ 0 46 4]]
image.png

Confusion matrixで分かるように青と緑はデータがかぶっているため誤判定がみられる。

ついでに**Mini Batch K-Means**も見てみる。
K-Meansを使うと計算コストが問題になるほどデータ数が多いときに代用するアルゴリズム
任意の数(mini batch)のSampleでK-Means同様にcentroidを移動させるが、mini batchを変えるたびに過去のcentroidと平均する。
精度はK-Meansより劣化するが無視できるほどらしい。

## Mini Batch K-Means
mb_kmeans = MiniBatchKMeans(n_clusters=3, random_state=0, batch_size=6).fit(x_train)
y_pred = kmeans.labels_
y_cent = kmeans.cluster_centers_

# 精度評価結果
print_evaluate(x_train, y_test, y_pred, 'macro')

colors_pred = [color_codes[x] for x in y_pred]
plt.scatter(x_train[:,0], x_train[:,1], c=colors_pred)
plt.scatter(y_cent[:,0], y_cent[:,1], s = 250, marker='*',c='y')
plt.show

accuracy = 0.37333333333333335
precision = 0.37286324786324787
recall = 0.37333333333333335
f1 score = 0.3730825663598773
confusion matrix =
[[50 0 0]
[ 0 2 48]
[ 0 46 4]]
image.png
もともとデータ数が少ないので精度が悪化はわからない。

3. GMMによるクラスタリング

では、まず**Gaussian Mixture Model (GMM)**でK-Meansのようにクラスタリングする。
最初にTrainingを行い、GMMはProbablistic ModelなのでPredictする。
もう少し言うと、predict()はラベルの予測を行い、predict_prob()はそのラベルとなる確率を出力する。
K-Meansはユークリッドな距離でクラスタ分けを行うので確率の概念がなかった。
GMMではこの確率を用いて後述のAnomly ditectionを行うこととなる。
まずモデルを作成する。

from sklearn.mixture import GaussianMixture

## Gaussian Mixture Model
gm = GaussianMixture(n_components=3, n_init=10)
gm.fit(x_train)

# fit()で収束したか確認
print(f'Convergence: {gm.converged_}')

y_pred = gm.predict(x_train)
y_pred = np.array(y_pred, dtype=np.float)  # Predict()の戻りの型がy_testと違うので一応合わせる。
y_pred_prob = gm.predict_proba(x_train)

# 精度評価結果
print_evaluate(x_train, y_test, y_pred, 'macro')

# 共分散行列
print('Covariance matrix = \n', gm.covariances_)

Convergence: True
accuracy = 0.02
precision = 0.019230769230769232
recall = 0.02
f1 score = 0.0196078431372549
confusion matrix =
[[ 0 50 0]
[ 1 0 49]
[47 0 3]]
ラベルの誤判定がK-Meansと比べ、2か所改善していることがわかる。

Covariance matrix =
[[[0.07200287 0.04402642]
[0.04402642 0.30034404]]
[[0.01088496 0.00594793]
[0.00594793 0.02955684]]
[[0.04585075 0.08489919]
[0.08489919 0.24667113]]]
n_components=3で指定したので3クラスタそれぞれの共分散行列(x_trainは2次元なので2行2列)を表示した。
おそらく、標準偏差、共分散ともに小さな真ん中の行列が、分布が比較的集中している緑のアヤメを示していると思う。

gm.score_samplesで等高線を追加した。

colors_pred = [color_codes[x] for x in y_pred]

# https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_pdf.html#sphx-glr-auto-examples-mixture-plot-gmm-pdf-py
x = np.linspace(0., 3.)
y = np.linspace(0., 7.)
X, Y = np.meshgrid(x, y)
XX = np.array([X.ravel(), Y.ravel()]).T
Z = -gm.score_samples(XX)
Z = Z.reshape(X.shape)

cont = plt.contour(X, Y, Z, norm=LogNorm(vmin=1.0, vmax=100.0),
                    levels=np.logspace(0, 3, 20), 
                    linestyles='dashed', linewidths=0.5)

cont.clabel(fmt='%1.1f', fontsize=12)
plt.xlabel('X', fontsize=12)
plt.ylabel('Y', fontsize=12)

plt.title('Negative log-likelihood predicted by a GMM')
plt.scatter(x_train[:,0], x_train[:,1], c=colors_pred)
plt.show

Screenshot from 2021-09-04 17-35-55.png

4. GMMによるクラスタ数のもとめ方

クラスタごとのGaussian disributionを現すパラメータは少なければ少ないほどよく、一方そのパラメータのlikelihood(尤度)は大きければ大きいほどよい。
とてもおおざっぱだけど、それを表現する基準がBIC(Bayesian Information Criterion)とAIC(Akaike Information Criterion)である。
詳細は上述Hands-on教科書のP267~269に解説あり。
クラスタ数を振ってBIC, AICをプロットすると最適なクラスタ数を得ることができる。

n_components = np.arange(1, 10)
models = [GaussianMixture(n, n_init=10).fit(x_train)
         for n in n_components]

fig, ax = plt.subplots(figsize=(9,7))
ax.plot(n_components, [m.bic(x_train) for m in models], label='BIC')
ax.plot(n_components, [m.aic(x_train) for m in models], label='AIC')

plt.legend(loc='best')
plt.xlabel('n_components')

image.png
ゑっ? よく見るチャートぢゃない・・・

サンプル数150が少ないのかもしれない。
GMMはProbabilistic modelなのでsample()メソッドで確率密度に沿ったデータの増産することができる。
そこで合計1000個(+850個)まで増やして再描画してみる。

x_new, y_new = gm.sample(850)
x_new = np.array(x_new).reshape(850,2)
y_new = np.array(y_new)
x_augment = np.vstack([x_train, x_new])
y_augment = np.append(y_test, y_new)

n_components = np.arange(1, 10)
models = [GaussianMixture(n, n_init=10).fit(x_augment)
         for n in n_components]

fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(n_components, [m.bic(x_new1) for m in models], label='BIC')
ax.plot(n_components, [m.aic(x_new1) for m in models], label='AIC')

plt.legend(loc='best')
plt.xlabel('n_components')

image.png
3の時がBIC,AICともに最小であることから、クラスタ数=3が尤もらしいとわかる。
(よかった)

もちろんK-Meansでもsilhouette_score()silhouette diagramで最適なクラスタ数を推定することができる。

5. GMMによるAnomaly Detection

調べたいデータがanomalyか否かは上述のGMMモデルのスコアとしてscore_samples()で簡単に調べることができる。
これで得られるスコアはデータの確率密度の濃さを示している。
なお、scikit-learnのscore_samples()は自然対数で返すためexp()した。

例えばデータ群の中の(0.2, 0.1)とデータ群から外した(2.5, 0.5)の2データについてスコアを見てみる。

x_compare = np.array([0.2, 1.0, 2.5, 0.5]).reshape(2, 2)
print(x_compare)

score_compare = gm.score_samples(x_compare)
print(np.exp(score_compare))

[[0.2 1. ]
[2.5 0.5]]
[7.56273077e-02 1.36493463e-24]

期待通り、前者が大きくなり、後者はほぼゼロとなっている。データ分布がGaussianなら色々使えそう。
もちろん、特異点の判断基準はドメインの要求によるため言及できない。

以上

23
23
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
23
23

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?