未経験からデータサイエンティスト

未経験からデータサイエンティストの勉強したことの備忘録とか雑記とか

機械学習による異常検知 2−1. 1変数正規分布

どうも、sue124です。

今回は前回に引き続き書籍「入門 機械学習による異常検知」のことを書いていこうと思います。
今回は1変数の正規分布の場合の異常検知について、異常度の閾値設定するまでの数式導出と、それを使ったPythonでの異常検知をやっていきたいと思います。

1.異常検知の閾値設定までの数式導出

以前の記事で紹介した通り、以下の3ステップで異常検知の閾値設定を行います。
sue124.hatenablog.com


ステップ1:分布推定
ステップ2:異常度の算出
ステップ3:閾値の設定

以下では正常とわかっているN個のデータ D = \{ x^{(1)}, x^{(2)}, ...,x^{(N)} \} から分布推定をして、新しく観測したデータ x' の判定をするケースを考えます。
それぞれのステップを詳しく見ていきます。

1ー1.分布推定

ここではデータの性質に応じた適切な確率分布のモデルを仮定します。
今回は確率分布のモデルに、平均値を中心として左右最小な分布となっている、正規分布を仮定して進めます。
1変数の場合の正規分布の式は、平均を \mu 、分散を \sigma^{2} とすると以下の通りです。

{ \displaystyle
N (x | \mu, \sigma^{2}) \equiv \frac{1}{(2\pi\sigma^{2})^{1/2}} \exp\left(-\frac{1}{2\sigma^{2}}(x - \mu)^{2}\right)
}

正規分布の式の導出は省略します。
(書籍の付録に、ラグランジュの未乗数法を使った導出が載っています。)

次に確率分布の式に出てくるパラメータ \mu\sigma を求めます。
これらのパラメータは通常、データに対して最も当てはまりが良い値を求める(最尤推定する)のですが、今回は1変数の場合ですので、一般に知られている平均値と標本分散の式から、以下の通り求めます。

{ \displaystyle
\hat{\mu} = \frac{1}{N} \sum_{n=1}^{N} x^{(n)}
}

{ \displaystyle
\hat{\sigma}^{2} = \frac{1}{N} \sum_{n=1}^{N} \left(x^{(n)} - \hat{\mu}\right)^{2}
}

ここで ^(ハット)はデータから推定した値であるということを示します。

1ー2.異常度の算出

新たな観測値 x' の異常度 a(x') を求めます。
それにあたってまず、1で求めた確率分布の負の対数尤度を取ると、以下の通りとなります。

{ \displaystyle
\frac{1}{2\hat{\sigma}^{2}} (x' - \hat{\mu})^{2} + \frac{1}{2}\ln({2\pi\hat{\sigma}^{2}})
}

x' に依存しない第二項を無視して、式の形を整えたものを以下の通り異常度と定義します。

{ \displaystyle
a(x') \equiv \left( \frac{x^{} - \hat{\mu}}{\hat{\sigma}} \right) ^{2}
}

上式の意味を考えると、以下の通りとなります。

  • 分子は標本平均からのずれ(=分布の中心からの距離)であり、これが大きくなるほど異常度が高い
  • 分母は標本標準偏差であり、元々のばらつきが大きい場合は多少のずれでは異常度は大きくなりにくく、ばらつきが小さいと、少しのずれで異常度が大きくなる

1ー3.閾値の設定

閾値を設定するにあたって、まずは異常度  a(x') の確率分布を考えます。
ここで下記の「ホテリング統計量の分布」の定理を使うと、 N \gg 1 の場合、異常度  a(x') は自由度1、スケール因子1のカイ二乗分布に従うことがわかります。

<定理:ホテリング統計量の分布>
1次元の観測データ D = \{ x^{(1)}, x^{(2)}, ...,x^{(N)} \} が独立に同じ分布 N(\mu, \sigma^{2} に従い、新たな観測値 x' も同じ分布に従うする。
このとき、異常度  a(x') の定数倍は自由度 (1, N-1) のF分布に従う。

{\displaystyle
\frac{N-1}{N+1} a(x') \sim F(1, N-1)
}

特に、 N \gg 1 の時は、  a(x') は自由度1、スケール因子1のカイ二乗分布に従う。

{\displaystyle
a(x') \sim \chi^{2}(1, 1)
}

書籍内では、異常度  a(x') の分母と分子がそれぞれカイ二乗分布に従って、その比がF分布に従うことの証明がありましたが、ここでは省略します。

カイ二乗分布の性質は既によく知られているので、それを活用します。
Rであれば、次のコードでカイ二乗分布のグラフを描画できます。

curve(dchisq(x, 1), xlim=c(0,4))

初めてRを使ってみたのですが、たったの1行でグラフを出力できることに驚きました。
出力結果は、下図のようになります。

f:id:sue124:20210131150347p:plain
カイ二乗分布

新たに観測した値の異常値  a(x') が異常かどうかを考えるには、それを正常時の分布に照らしてどの程度「あり得ない」のかを計算する必要があります。
そのためには、「その異常値より右側の面積が全体の面積( a(x') の全区間積分)に占める割合」を求めます。
 a(x') = 2 の場合ですと、下図の赤斜線部の面積の割合を求めることになります。

f:id:sue124:20210131151950p:plain
 a(x') = 2 の場合

これによって、物理量などの特性に依存せずに客観的に異常判定の閾値を決めることができます。

ここで話を異常度の閾値の設定に戻します。
例えば、「閾値を1%値に選ぶ」場合、「赤斜線部の面積が全体の1%になるような  a(x')閾値にする」ということになります。

これで異常度がこの値より大きい場合は、それを「異常」と判定することができるようになりました。

2.Pythonでの異常検知

では上記を活かしてPythonで異常検知をしていきたいと思います。
データセットとしては、irisを使っていきます。
まずはデータを準備します。

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from scipy.stats import chi2

iris = sns.load_dataset('iris')

# 種ごとにDataFrameを分割する
setosa = iris[iris['species']=='setosa']
versicolor = iris[iris['species']=='versicolor']
virginica = iris[iris['species']=='virginica']
# setosaのデータに他の種のデータを1つずつ混入させる
df = pd.concat([setosa, versicolor.head(1), virginica.head(1)]).reset_index()
display(df.tail())

setosaのデータに他の種のデータを1つずつ混入させたDataFrameの末尾はこんな感じになっています。

f:id:sue124:20210211230144p:plain

ヒストグラムを見てみます。

plt.hist(df['sepal_length'])
plt.xlabel('sepal length (cm)')
plt.ylabel('count')
plt.show()

混入させたデータは、元のヒストグラムの山から外れた位置にいます。
※setosaデータの山が左右対象でないので厳密には正規分布を当てはめるのは本当は相応しくないのですが、正規分布以外の場合はまた別の機会に書きます。

f:id:sue124:20210211230442p:plain


データを分布推定用のものと、異常度評価用に分けて異常度を算出していきます。
異常度の閾値は1%値にすることにします。

 データ分割
df_train = df[:45]
df_test = df[45:].reset_index(drop=True)
# %%
# 異常度算出
mu = df_train['sepal_length'].mean()
s2 = ((df_train['sepal_length'] - mu)**2).mean()
a = (df_test['sepal_length'] - mu)**2 / s2

print(f'標本平均:{mu:.3g}')
print(f'標本分散:{s2:.3g}')
print(f'異常度:\n{a}')

# 閾値の設定
th = chi2.isf(0.01, 1)
print(f'カイ2乗分布による1%水準の閾値: {th:.3g}')
# %%
plt.plot(range(len(a)), a, 'bo')
plt.axhline(y=th, color='red')
plt.ylabel('anomaly score')
plt.show()

出力はこのようになり、最後に付け足した種が異なる2つがうまく異常と判別できていることがわかります。
このようにすることで異常なデータを見つけることができるので、測定ミスなどで混入した不要なデータを取り除く作業(クレンジング)に使うことができます。

標本平均:5.01
標本分散:0.129
異常度:
0     0.346715
1     0.061468
2     1.314829
3     0.649251
4     0.000960
5    30.773146
6    12.923550
カイ2乗分布による1%水準の閾値: 6.63

f:id:sue124:20210211231411p:plain

3.まとめ

今回は1変数の正規分布の場合の異常検知についてを書きました。
数式を導出し、Pythonでirisデータセットを使って、異常検知を順を追ってみていきました。
次回はこれを多変数の場合に拡張したいと思います。
お読みいただきありがとうございました。