【Pythonで異常検知】Chapter 1. 1変数正規分布に基づく異常検知

Python
Sponsored

概要

この章では、以下の手順にしたがって1変数データの異常検知をPythonで実践することを目標とする。

  1. 訓練用データを正規分布にフィッティングする
  2. 得られた正規分布からテスト用データの異常度を求める
  3. 異常度の閾値を設定し、それを上回るデータを異常と判定する

データの準備

異常検知を実践するために、今回はScikit-learnのbrest_cancerのデータセットを用いる。このデータセットは、乳房にできた腫瘍の様々な特徴量の値と、その腫瘍が悪性・良性のいずれであったかを記載したラベルから成る。

変数dataにデータセットを読み込むと、data["data"]に特徴量の値、data["target"]に特徴量と関連付けられたラベルが格納される。ラベルは0が悪性、1が良性を意味し、特徴量の値からこれらのラベルを予測することが課題となる。

import numpy as np
from sklearn import datasets

data = datasets.load_breast_cancer() # 乳癌データ読み込み
target_0 = data["target"][data["target"] == 0] # 悪性のラベルのみ集める(要素はすべて0)
target_1 = data["target"][data["target"] == 1] # 良性のラベルのみ集める(要素はすべて1)

このデータセットは30の特徴量を持つため、今回はその中からworst areaの特徴量だけを用いる。また、このデータセットを訓練用とテスト用に分割することになるが、訓練用には良性(正常)のデータのみを用いるため、良性212例を訓練用としてモデルを作成し、残りの良性145例と悪性212例に対するラベルの予測を行うこととする。

data_w = data["data"][:,np.where(data["feature_names"] == "worst area")[0][0]] # worst areaのデータのみ抜き出す
data_w_0 = data_w[data["target"] == 0] # 悪性のworst areaデータ(212例)
data_w_1 = data_w[data["target"] == 1] # 良性のworst areaデータ(357例)

X_train = data_w_1[:212] # 訓練用に212例の良性データを用いる
X_test = np.hstack((data_w_0, data_w_1[212:])) # テスト用に212例の悪性データ+145例の良性データを用いる
y_test = np.hstack((target_0, target_1[212:])) # テストデータに付けられるラベル

正規分布へのフィッティング

訓練用データの特徴量のばらつきを、ヒストグラムで表示すると以下のようになる。

from matplotlib import pyplot as plt
plt.hist(X_train, bins=20) # 20個区切りでworst areaのヒストグラムを作成
plt.title("Histgram")
plt.xlabel("worst area")
plt.ylabel("Number")
plt.show()

特徴量の分布は500付近で多く、そこから遠ざかるにつれて数が少なくなる山型の分布を取っている。今回の異常検知では、この特徴量の確率分布が正規分布にしたがうと考えて、この山型が最も適合する正規分布の形へフィッティングする。このとき、得られるのはフィッティングした正規分布の平均\(\mu\)と標準偏差\(\sigma\)である。なお、正規分布の導出とその性質についての解説は以下の記事を参照のこと。

ここでは、Scipyに実装されているnorm.fitを用いてフィッティングを行う。

from scipy.stats import norm
param = norm.fit(X_train) # mu=548.7400943396226, sigma=158.27941541866434

x = np.linspace(200,1200,100) # x軸の値を作成
pdf_fitted = norm.pdf(x, loc=param[0], scale=param[1]) # paramから作成した正規分布の確率密度関数
plt.hist(X_train, bins=20, density=1) # ヒストグラムを確率密度関数化して表示
plt.plot(x,pdf_fitted) # 確率密度関数を表示
plt.title("Histgram & Norm")
plt.xlabel("worst area")
plt.ylabel("Probability Density Function (PDF)")
plt.show()

上図のオレンジ色の曲線が得られた正規分布である。なお、上図では特徴量の分布も確率密度として表示し直している。フィッティングの結果\(\mu\simeq 549, \sigma\simeq 158\)が得られ、このことから、特徴量の値が500程度になることはよくあることだが、100や1,000になることはほとんどあり得ないことが読み取れる。

異常度の計算

この、特徴量の値の起こりやすさを数値的に表現したものが「異常度」である。1変数正規分布では、特徴量が値\(x\)を取ることについての異常度\(a(x)\)は次のように表される。

この値は\(x\)が\(\mu\)に近いときに小さく、\(\mu\)から遠いときに大きい。すなわち、名前の通り異常さの度合いを表している。また、データ数(今回は症例数)が非常に大きい時は、この異常度\(a(x)\)自身が自由度1、スケール因子1の\(\chi^2\)分布に従う。

この\(\chi^2\)分布を表示したのが下図である。

from scipy.stats import chi2
x = np.linspace(0,5,100) # x軸の値を作成
chi_2 = chi2.pdf(x, 1, scale=1) # 自由度1, スケール因子1のχ2乗分布
plt.plot(x,chi_2) # χ2乗分布の表示
plt.title("Chisquare")
plt.xlabel("Abnomarity")
plt.ylabel("Probability Density Function (PDF)")
plt.show()

閾値の設定

上図に表示した\(\chi^2\)分布は確率密度関数であるため、これを\((0,\infty)\)の区間で積分すると1になる。また、図によると異常度が0~4程度の範囲が大部分の面積を占めているため、異常度が4以上となる確率はほとんどないことがわかる。異常検知の最終ステップでは、「異常度が〇以上になることは、△%程度の確率でしか起こり得ない(すなわち、ほとんど起こり得ない=異常である)」という異常度の閾値を設定し、それを超える異常度を持つデータ(症例)を異常(悪性)であると判断する。

たとえば、起こり得る確率が5%未満となる異常度の値は次のように求められる。

th_5 = chi2.ppf(0.95, 1, scale=1) # 3.8414588206941236

x_th_5 = x[x >= th_5] # x軸で閾値を超える範囲
y0 = np.zeros_like(x_th_5) # y=0
plt.plot(x, chi_2) # χ2乗分布の表示
plt.fill_between(x_th_5, y0, chi_2[-1*len(y0):], where=y0<chi_2[-1*len(y0):], color="r") # "異常"と判定される領域を塗りつぶす
plt.title("Chisquare")
plt.xlabel("Abnormality")
plt.ylabel("Probability Density Function (PDF)")
plt.show()

\(th_5\simeq 3.84\)より、異常度が3.84を超える確率は5%程度である。正常例の5%程度しか取り得ない値の特徴量であればそれはもう異常であると判断してもよく、ここでは上図の赤色で塗りつぶした範囲の異常度をもつデータはすべて異常とするモデルを作成する。

def predict(X, param, th):
    a = ((X - param[0])/param[1])**2 # 異常度の計算
    pred = np.zeros_like(a) # 予測したラベル(すべて0(悪性)として初期化)
    pred[a <= th] = 1 # 異常値が閾値を下回る場合は、ラベルを1(良性)とする
    return pred # 予測を返す

from sklearn.metrics import confusion_matrix
print(confusion_matrix(y_test, predict(X_test, param, th_5)))

# 177  35
#   7 138
 予測結果
悪性良性
真の状態悪性17735
良性7138

predict関数はデータ、パラメータ(\(\mu\)と\(\sigma\))、閾値を引数に与えることで、ラベルの予測結果を返す。また、sklearn.metrics.confusion_matrixは予測結果と真の状態の混同行列を返す。

今回は悪性例212例のうち177例が正しく悪性と判断され、35例の悪性は見落とされた。

ここで、閾値を変えることによって結果を改善することができる場合がある。次の実行例では、異常度の閾値を10%未満と、異常の範囲を広く取っている。

th_10 = chi2.ppf(0.90, 1, scale=1) # 2.705543454095417
print(confusion_matrix(y_test, predict(X_test, param, th_10)))

# 186  26
#  19 126
 予測結果
悪性良性
真の状態悪性18626
良性19126

その結果、正しく悪性と判定される割合が増加したが、今度は先程よりも多くの良性例が間違って悪性と判断されていることがわかる。

次は逆に、異常度の閾値を1%未満とし、異常の範囲を狭めてみる。

th_1 = chi2.ppf(0.99, 1, scale=1) # 6.634896601021213
print(confusion_matrix(y_test, predict(X_test, param, th_1)))

# 161  51
#   1 144
 予測結果
悪性良性
真の状態悪性16151
良性1144

今度は間違って悪性と判定される割合が非常に少なくなったが、かなり多くの悪性例を見落としている。

閾値をどのように定めるかによって予測性能は大きく変化しうるが、どの閾値が最適であるかの評価は難しい。最終的には異常検知を行う理由こそが閾値を決定するものとなり、たとえば今回の乳癌の判定のように異常(悪性)を見落とすと大きな不利益を被る場合には、多少偽陽性が増加しても異常の範囲は広く取った方がよいと考えられる。

Comments