異常検知の基礎
2018/01/11
  • このエントリーをはてなブックマークに追加

異常検知の基礎

はじめに

カブクで機械学習エンジニアをしている大串正矢です。今回は異常検知の基礎について書きます。

背景

異常検知の手法は多種多様に存在していますがウェブ上にまとまった情報が日本語でないため記述することにしました。ただ全ての内容をこの記事で記述すると長くなるため今回は基礎的な内容にフォーカスして記述します。この記事の内容をベースに他の異常検知に関する記事を記述する予定です。

異常検知の基本的なアプローチ

下記の3つのフェーズで構成されています。
  • 分布推定
    • モデルを定義し、正常なデータから学習します。
  • 異常度の定義
    • 上記のモデルからのずれの度合いである異常度を定義します。
  • 閾値の設定
    • 異常度がある値より大きければ異常と判定できるような閾値を設定します。
ここからは簡単なケースとして正規分布を想定して書いていきます。

分布推定

正規分布の数式は下記のようになります。
\begin{align}
N(x|\mu, \sigma^2) = \frac{1}{(2\pi\sigma^2)^{1/2}}\exp{\frac{-1}{2\sigma^2}(x-\mu)^2}
\end{align}

図で表すと下記のようになります。


上式のμ(平均)とσ(分散)が決まれば正規分布を求めることができます。
これは単純に下記の式で導出できます。

\begin{align}
\mu = \frac{1}{N}\sum^{N}_{n=1} x^n
\end{align}

\begin{align}
\sigma = \frac{1}{N}\sum^{N}_{n=1} (x^n - \mu)^2
\end{align}

上記によって導出したμ(平均)とσ(分散)によって正規分布の分布推定が終わります。
正規分布をモデルと定義し、平均と分散をデータから導出して適用しているので分布推定に当てはまります。

異常度の定義

負の対数尤度を異常度として採用します。平均と分散の導出に使用していない新たな観測値であるx’を使用して分布推定で使用した正規分布の対数を計算すると

\begin{align}
\frac{1}{2\sigma^2}(x' - \mu)^2 + \frac{1}{2}\ln(2\pi\sigma^2)
\end{align}

第2項は新たな観測値x’に依存しないので無視して、式の形を綺麗にするために2をかけて異常度a(x’)を求めます。
\begin{align}
a(x') = \frac{1}{\sigma^2}(x' - \mu)^2 = (\frac{x' - \mu}{\sigma})^2
\end{align}

上記で異常度の導出は終わりです。上式では平均が中心である正規分布において、その平均からの離れ度合いを標準偏差で割ることによって正規化しています。イメージとしてはデータに応じた中心からの離れ度合いになります。


閾値の設定


閾値の設定は簡単な手法として分位点の設定があります。
分位点を使用する場合は先ほど導出した異常度を全観測データに対して適用します。分位点が2%の場合はその中の2%は滅多に観測されないので異常とするという考え方です。
仮に全観測データが200件の場合はその全観測データに対して異常度を導出し、異常度の降順に並べ直します。
200 * 0.02=4なのでその中の上位4番目の異常度を閾値とする考え方です。

この手法自体はシンプルですがデータのばらつきに弱いです。
そこで異常度の確率分布から閾値を設定するための手法を紹介します。
異常度がどのような分布に従うか知る必要があります。ホテリング理論の定理2.1によると

1次元の観測データ

\begin{align}
D = {x^{(1)},....,x^{(N)}}
\end{align}

の各観測値が独立に同じ正規分布に従い、新たな観測値x'も同じ正規分布に独立に従うとする。このとき異常度a(x')の定数倍は、自由度(1, N-1)のF分布に従う。すなわち


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

特に N >> 1の時はa(x')そのものが自由度1、スケール因子1のカイ2乗分布に従う

\begin{align}
a(x') \sim \chi^2(1, 1)
\end{align}

この定理を用いて異常度が自由度1、スケール因子1のカイ2乗分布に従うと考え閾値設定をしていきます。

カイ2乗分布は実は正規分布との関係があります。実際に下記のコードを動作させてその関係を把握してみます。


import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
from matplotlib import animation as ani


k = 1
n = 30000
cum = np.zeros(n)
for i in range(k):
    # 標準正規分布に従う乱数を30,000個生成
    x = np.random.normal(0, 1, n)
    # 生成した乱数を2乗する
    x2 = x**2
    # 足し合わせた数が自由度となります。
    cum += x2
# ヒストグラムの描画
plt.figure(figsize=(7,5))
plt.title("chi2 distribution.[k=1]")
plt.hist(cum, 80, color="lightgreen", normed=True)
# 自由度1のカイ二乗分布の描画
xx = np.linspace(0, 25 ,1000)
plt.plot(xx, chi2.pdf(xx, df=k, scale=1), linewidth=2, color="b")
plt.show()


自由度1のカイ2乗分布は正規分布の2乗とほぼ等価であることが分かります。数学的な証明は入門-機械学習による異常検知―Rによる実践ガイド-井手-剛のp.30から記述されているので気になる方はそちらをご参照ください。
ではカイ2乗分布と異常スコアの関係に戻りましょう。
\begin{align}
a(x') \approx \chi^2(1, 1)
\end{align}
カイ2乗分布を数式で表すと下記になります。
\begin{align}
u = a(x')
\end{align}
\begin{align}
\chi^2(u|k, s) = \frac{1}{2s\Gamma(k/2)}\frac{\mu}{2s}^{(k/2-1)}\exp(-\frac{1}{2s})
\end{align}

kは自由度で分布の形状を変更する作用があります。観測する変数の数だけ自由度が増えます。sはスケール因子で分布の大きさを調整します。sが分母にあるため小さいほど分布が大きくなり、大きいほど小さくなります。これによりカイ2乗分布の面積を調整することができます。

Γはガンマ関数を表しています。
ここで重要なのがカイ2乗分布の面積が確率になることです。
異常度の変数は1つなので自由度は1でスケール調整の必要もないので1として面積を導出すると
\begin{align}
\alpha = \int_{ath}^{\infty} \chi^2(u|1, 1) du = 1 - \int_0^{ath} \chi^2(u|1, 1) du
\end{align}

αをここでは0.01にすると下記の図の面積が0.01になるような異常スコアの閾値athを導出することになります。
図で表すと下記のようになります。


Pythonによる実装


ここからはPythonとサンプルのデータを用いて具体的な実装に入っていきます。ガンマ分布はライブラリが用意されているため上記のような数式を記述せずとも使えます。
使用したデータセットはDavis.csvです。このデータの中の体重の項目を使用しているため明らかに体重が大きいもしくは小さいデータは異常と判定します。


下記がコードになります。


import csv
from scipy import stats
import argparse


def main():
   parser = argparse.ArgumentParser(description="hotelling theory")
   parser.add_argument("-d", "--data", metavar="data",
                       type=str, default='../data/raw/davis.csv',
                       help="setting test data")
   parser.add_argument("-t", "--threshold_rate", metavar="threshold_rate",
                       type=float, default=0.01,
                       help="setting threshold_rate")
   # parse arguments
   args = parser.parse_args()
   # データセットの読み込み
   num = []
   data = []
   with open(args.data, 'r', encoding="utf-8") as f:
       reader = csv.reader(f)
       header = next(reader)
       for row in reader:
           num.append(int(row[0]))  #標本番号を取得
           data.append(int(row[2])) #体重データを取得
   # 標本平均
   mean = np.mean(data)
   # 標本分散
   variance = np.var(data)
   # 異常度
   anomaly_scores = []
   anomaly_scores_dict = {}
   for x in data:
       anomaly_score = (x - mean)**2 / variance
       anomaly_scores.append(anomaly_score)
       anomaly_scores_dict.update({anomaly_score: x})
   # カイ二乗分布による1%水準の閾値
   threshold = stats.chi2.interval(0.99, 1)[1]
   for k, v in anomaly_scores_dict.items():
       if k > threshold:
           print("anomaly weight {0} kg, anomaly score {1}".format(anomaly_scores_dict[k], k))
   # 結果の描画
   plt.plot(num, anomaly_scores, "o", color = "b")
   plt.plot([0,200],[threshold, threshold], 'k-', color = "r", ls = "dashed")
   plt.xlabel("Sample number")
   plt.ylabel("Anomaly score")
   plt.ylim([0,100])
   plt.show()


if __name__ == '__main__':
   main()

下記が動作結果です。
ターミナルに異常な体重の値が出ます。
図は異常値のスコアと閾値です。



anomaly weight 119 kg, anomaly score 12.483415666901907
anomaly weight 166 kg, anomaly score 44.28387438249824


下記が動作結果です。青が異常スコアのデータで赤が閾値になります。先ほど異常と判定された166kgと119kgのデータのみ閾値を超えています。


今後

次回はこの記事を踏まえた上で各異常検知の手法のユースケースと特に時系列の異常検知に着目し、その内容について紹介します。

最後に

弊社では異常検知以外にも物体検出、3次元データ検索エンジンの開発をしています。これらに興味があるエンジニアがいらっしゃれば絶賛採用中なので是非、弊社へ応募してください。

参考


入門-機械学習による異常検知―Rによる実践ガイド-井手-剛
Pythonでお手軽・異常検知 [ホテリング理論編]
【統計学】正規分布とカイ二乗分布の関係を可視化してみる。