Naoism's Blog

データサイエンス、機械学習、競技プログラミング、音楽、勉強

Pythonで異常検知(1):1変数データに対して

お久しぶりです。

新年度が始まり、異常検知という新たな分析案件に取り組むことになりまして。
案件を進めていく前に、分野の基礎知識を学び、分析方針を立てることが大事たと考えているため、まず以下の本をもとに勉強を進めていくことを決めました。

www.coronasha.co.jp

「入門 機械学習による異常検知」です。
(MLPの方と迷ったのですが、井手先生のHPに上記の本の続編と位置づけられると記載してあったのでやめました。)

この本を買うときに迷ったのが、記載コードがRで書かれているということ。
個人的にpythonのほうが圧倒的に書き慣れているので、購入する際に少し躊躇してしまいました。

ですが、いざ買って読んでみるとメインは理論であり、それを補足する程度にRのコードが書かれているという形式でした。
また、その理論も丁寧で、買ってよかったと心から思える非常に素晴らしい本でした。

しかし、私と同じようにpythonのほうが好きで購入に踏み込めない、という方も多くいると思うので、勉強がてら本に記載されているコードをpythonで書いていこうと思います。

使用データ

今回は、井手先生本でも利用していた、Davis Datasetを用いていきます。
このデータには、性別、体重(実測)、身長(実測)、体重(自己申告)、身長(自己申告)という5つのデータが200人分記録されています。

データは下記からダウンロードしました。
https://vincentarelbundock.github.io/Rdatasets/datasets.html

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

data = pd.read_csv("Davis.csv", index_col=0)
print(data.head())

  sex  weight  height  repwt  repht
1   M      77     182   77.0  180.0
2   F      58     161   51.0  159.0
3   F      53     161   54.0  158.0
4   M      68     177   70.0  175.0
5   F      59     157   59.0  155.0

今回は、本データにおける1変数「体重(weight)」に対して異常検知分析をしていきます。
とりあえずどんな分布かを見るためにヒストグラムを書いてみましょう。

plt.figure(figsize=(12,8))
plt.hist(data["weight"], bins=20)
plt.xlabel("Weight (kg)", size=15)
plt.ylabel("Frequency", size=15)

f:id:mashkun:20200403231532p:plain

分布は、若干左右非対称ですが、概ねひと山の形となっています。

この山を表現できる確率分布はたくさんありそうですが、まずは正規分布を当てはめることを考えていきます。

ホデリング理論

ホデリング理論とは、観測データが正規分布に従うことを仮定した下で、分布から逸脱する異常値を検出するという方法論です。

以下計算手順は井出先生本の引用です。詳細を知りたい方はぜひ買ってみてください。


 手順 (ホデリングT^ 2法(1次元))
 (0) 準備
  異常が含まれていないか、含まれていたとしてもごく少数と思われるデータセットを用意する。
  異常判定の閾値を確率値αで与え(例えば0.01や0.03)、
  カイ二乗分布の表から、異常度の閾値α _ {th}を求めておく。

 (1) ステップ1(分布推定)
  標本平均および標本分散を計算する。

 (2) ステップ2(異常度の計算)
  新たな観測値x'が得られるたび、異常度
  a(x')=(\frac{x′−u}{\hat{σ}})^ 2
  の値を計算する。

 (3) ステップ3(閾値判定)
  異常度が閾値α _ {th}を超えたら異常と判定する。


本では、Rを用いて上記手順を試していますが、今回はpythonでやってみます。

手順(0)

データセットは上記のDavisを用います。
また、異常判定の閾値は1%とします。
異常度の閾値は、カイ二乗分布を表から求めてもいいですが、せっかくpythonを使うので、stats.chi2.interval()を使って計算することにしましょう。

手順(1)

単純に標本平均、標本分散を求めます。

mu = np.mean(data["weight"]) # 標本平均
s2 = np.var(data["weight"]) # 標本平均
print("標本平均:", mu)
print("標本分散:", s2)

標本平均: 65.8
標本分散: 226.71999999999977
手順(2)

実際は新たな観測値に対して、異常度を計算しますが、
ここでは簡単のため、分布推定に用いた元データの各観測値について異常度を計算していきます。
(ホテリング理論の前提である独立性の過程を破っているので、厳密に言えば誤用。)

a = (data["weight"] - mu)**2 / s2 # 異常度
th_01, th_99 = stats.chi2.interval(0.98, 1) # (alpha, k) 自由度kのカイ2乗分布でalpha×100%となる範囲を取得
手順(3)

上記手順で求めた異常度と閾値を可視化してみます。
その下に、どのサンプルが異常と判断されたかを確認するために、元データのグラフも載せておきます。

plt.figure(figsize=(12,8))
plt.scatter(x=np.arange(len(a)), y=a)
plt.plot([0,200],[th_99, th_99], color = "r", ls = "dashed")
plt.xlabel("Sample number", size=15)
plt.ylabel("Anomaly score", size=15)
plt.show()

plt.figure(figsize=(12,8))
plt.plot(data["weight"])
plt.xlabel("Sample number", size=15)
plt.ylabel("Weight (kg)", size=15)
plt.show()


f:id:mashkun:20200403231708p:plain

f:id:mashkun:20200403232101p:plain


図より、異常と判断されたのは2点であり、元データの中で160kgと明らかな外れ値の異常度がより大きく計算されていることがわかります。

ここまでがホデリング理論に沿った異常検知の方法です。

本当に正規分布で仮定していいの?

今までデータを正規分布で仮定して分析をしていきました。
しかし、上記データは左右非対称であり、体重とは非負の値です。
それにも関わらず、正規分布を仮定しても良いのでしょうか?

実はこのような、①正の値しか取らない、②分布が左右対称ではない、という特徴を持つデータに対しては、カイ二乗分布、またはガンマ分布の当てはめが有効と本の3章に書かれています。

ガンマ分布を当てはめて分析

では、同じデータにガンマ分布を当てはめてみましょう。

ガンマ分布は以下の式で表されます。

G(x|k,s)\equiv\frac{1}{s\Gamma(k)}(\frac{x}s)^ {k-1}exp(-\frac{x}s)


ガンマ分布のパラメータを算出する必要がありますが、ここではモーメント法(と最尤法)により計算します。
計算の詳細を学びたい方は是非本をご購入ください。

以下ガンマ分布による異常検知の手順です。


 手順 (ガンマ分布による異常検知)

 (1) ステップ1(分布推定)

  モーメント法、あるいは数値最適化による最尤法によって、ガンマ分布のパラメターを推定する。

 (2) ステップ2(異常度の計算)
  新たな観測値x'を得るたび、異常度
  a(x')=\frac{x'}{\hat{s}}-(\hat{k}-1)\ln{\frac{x'}{\hat{s}}}
  を計算する。

 (3) ステップ3(閾値判定)
  訓練データを使ってあらかじめαパーセンタイルに対応する異常度の値α _ {th}をを求めておき、
  a(x')>α _ {th}(a)なら警報発報。

手順(1)

まず準備のため、標本数、標本平均、標準偏差を求めます。

N = len(data["weight"]) # 標本数
mu = np.mean(data["weight"]) # 標本平均
si = np.std(data["weight"]) # 標準偏差
print("標本数:", N)
print("標本平均:", mu)
print("標準偏差:", si)

標本数: 200
標本平均: 65.8
標準偏差: 15.057224179774963

その後、モーメント法を用いてksの値を計算します。
最尤法の計算結果とも比較したかったので念の為計算してみました。

kmo = (mu/si)**2  # モーメント法によるkの推定値
smo = si**2/mu  # モーメント法によるsの推定値
kml, loc, sml = stats.gamma.fit(data["weight"], floc=0)  # 最尤法によるkとsの推定値

print("モーメント法によるkの推定値:", kmo)
print("モーメント法によるsの推定値:", smo)
print("最尤法によるkの推定値:", kml)
print("最尤法によるsの推定値:", sml)

モーメント法によるkの推定値: 19.09685956245591
モーメント法によるsの推定値: 3.4455927051671704
最尤法によるkの推定値: 22.48530359200564
最尤法によるsの推定値: 2.926355863097812
手順(2)

計算式に従って、異常度を計算していきます。

a = (data["weight"] / smo) - (kmo - 1)*np.log(data["weight"] / smo)  # 異常度
手順(3)

ここでは、上位1%に当たる分位点を閾値としています。
図の点線がその閾値です。

th = sorted(a, reverse=True)[int(0.01*N)-1]

plt.figure(figsize=(12,8))
plt.scatter(x=np.arange(len(a)), y=a)
plt.plot([0,200],[th, th], color = "r", ls = "dashed")
plt.xlabel("Sample number", size=15)
plt.ylabel("Anomaly score", size=15)

f:id:mashkun:20200404184700p:plain


図を見るとわかるように、異常度の値が正規分布のときと異なるので気をつけましょう。
これは、仮定した分布によって異常度の算出式が異なるためです。

正規分布の結果と比較してみると、ガンマ分布で当てはめたときのほうが160kgの人の異常度が相対的に小さくなっていることがわかります。
これは、正規分布が左右対称のデータを想定しているのに対し、ガンマ分布が右に偏った分布を許しているからだと思われます。

右に大きく歪んでいるデータもまぁあり得るよね。というガンマ分布の気持ちが汲み取れますね。

最後に

今回は、1変数データに対して正規分布、ガンマ分布を仮定して異常診断をしてみました。

異常検知については実務で使うので日々勉強を続け、学んだことをこのブログに書いていこうと思います。
初めて勉強する分野なので間違ったことを記載してしまうかもしれませんがそこはご了承を。。