1次元正規分布に基づく異常検知の理論とPythonによる実装

はじめに

　異常検知とは、大多数のデータとは振る舞いが異なるデータを検出する技術のことです。異常検知は、膨大なデータが収集可能となった現代におけるデータ活用のひとつとして脚光を浴びています。

統計的異常検知の考え方

　異常検知にもいろいろアプローチがあります。例えば、距離に基づくアプローチ、統計学に基づくアプローチ、機械学習に基づくアプローチなどが挙げられます。今回はその中でも統計学に基づくアプローチについて書きます。統計学に基づく異常検知のアプローチは統計的異常検知とも呼ばれます。
　統計的異常検知とは、データの生成機構が確率モデル（確率分布）で表現できると仮定した場合の異常検知の方法論です。これはデータから得られる「知識」を確率モデルという形で表現していると捉えることもできます。

統計的異常検知の手順

統計的異常検知は基本的に以下の3ステップからなります。

観測データからデータ生成の確率モデルを学習
さらに細かく分けると以下の２ステップあります。
① 未知パラメータを含む確率分布を仮定
② データから未知パラメータを推定
学習したモデルを基に、データの異常度合いをスコアリング（異常度の定義）
閾値の設定

今回は、統計的異常検知の中でも最も基本的な、1次元正規分布に基づく異常検知について、上記の３ステップを一つずつ解説していきます。

データの準備

用いるデータは、Rのデータセットにある「Davisデータ」というもので、200人の性別、測定した身長・体重、および自己申告の身長・体重が記録されています。 https://vincentarelbundock.github.io/Rdatasets/datasets.html

　前提として、今回扱うデータは測定した体重だけであり、１変数です。すなわち確率的に揺らぐ値（確率変数）が1種類だけ存在するということです。もちろん、これから示す方法論は、多変量にも拡張できますが、1変数から始めた方が数式が追いやすいです。
　横軸を標本番号、縦軸を体重として、データを可視化したものが以下になります。

f:id:hirotsuru314:20170919224410p:plain:w550

　また、このデータのヒストグラムは下のようになります。

f:id:hirotsuru314:20170919224425p:plain:w550

　どちらもパッと見で異常値が見つけられると思います。このデータの異常判定をしたいとき、「120kg以上は異常だ！」としてif文で検知するのは簡単ですが、その120kgという値には何の根拠もなく、説得力に欠けると思います。これを統計学を駆使して、より客観的な水準で判定してやろうというのが今回の話です。
観測データがN個あるとき、データをまとめて $D$ という記号で表すとします。

$\displaystyle D= \{x_1,x_2,...,x_N\}$

ここで、観測データが多次元の場合は、 $x$ が列ベクトルとして表されます。また、このデータの中には異常な観測データが含まれていないか、含まれていたとしてもその影響は無視できると仮定します。
これから前述した統計的異常検知の3ステップに沿って話を進めていきます。

1. 観測データからデータ生成の確率モデルを学習する

① 未知パラメータを含む確率分布を仮定
　上記のヒストグラムによると、体重データの分布は若干左右非対称ではありますが、おおむね山の形になっています。そこで、それぞれの観測データ $x$ が平均 $\mu$ 、分散 $\sigma^2$ の正規分布に従うと仮定します。

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

② データから未知パラメータを推定
　正規分布に含まれているパラメータは平均 $\mu$ 、分散 $\sigma^2$ の二つです。言い換えると、この二つのパラメータが決定すると、分布の形が一意に決まるということです。これらのパラメータを観測データから推定します。パラメータ推定には最尤推定を用います。最尤推定は、「観測データが得られる確率」をパラメータの関数とみなした尤度関数を最大化するようにパラメータを決定する推定法です。

$\displaystyle \begin{eqnarray} p(D|\mu, \sigma^2) & = & \prod_{n=1}^N N(x_n|\mu, \sigma^2) \\ & = & \prod_{n=1}^N\frac{1}{\sqrt{2\pi \sigma^2}}\exp{\left\{-\frac{(x_n-\mu)^2}{2\sigma^2}\right\}} \\ & = & (2\pi \sigma^2)^{-\frac{N}{2}}\exp{\left\{-\frac{1}{2\sigma^2}\sum_{n=1}^{N}(x_n-\mu)^2 \right\}} \end{eqnarray}$

ここで、計算の都合上、尤度関数の自然対数をとった対数尤度関数を最大化します。対数は単調増加であるため、尤度関数の最大化は対数尤度関数の最大化と同値です。
$\beta=\frac{1}{\sigma^2}$ とした場合、対数尤度尤度関数は、

$\displaystyle \ln p(D|\mu, \sigma^2)=\frac{N}{2}\ln \beta - \frac{N}{2}\ln 2\pi - \frac{\beta}{2}\sum_{n=1}^{N}(x_n-\mu)^2$

これを最大化するパラメータを求めるため、 $\mu$ と $\beta$ でそれぞれ偏微分してゼロと等置します。

$\displaystyle \frac{\partial}{\partial \mu}\ln p(D|\mu, \sigma^2)=0 \tag{1}$

$\displaystyle \frac{\partial}{\partial \beta}\ln p(D|\mu, \sigma^2)=0 \tag{2}$

$(1)$ より

$\displaystyle \begin{eqnarray} \frac{\partial}{\partial \mu}\ln p(D|\mu, \sigma^2) & = & \beta\sum_{n=1}^{N}(x_n-\mu) \\ & = & \beta \biggl(\sum_{n=1}^{N}x_n-N\mu\biggl) \\ \end{eqnarray}$

$\displaystyle ∴　\mu=\frac{1}{N}\sum_{n=1}^{N}x_n$

$(2)$ より

$\displaystyle \frac{\partial}{\partial \beta}\ln p(D|\mu, \sigma^2) = \frac{N}{2\beta}-\frac{1}{2}\sum_{n=1}^{N}(x_n-\mu)^2$

$\displaystyle ∴　\sigma^2=\frac{1}{N}\sum_{n=1}^{N}(x_n-\mu)^2$

　平均・分散の推定値は私たちがよく知っている平均・分散の定義そのものだということがわかると思います。これらの結果は、平均・分散の推定値として、観測データの標本平均・標本分散を採用することを意味しています。また、推定した平均・分散にはデータからの推定値であることを明示するため「^（ハット）」をつけます。
したがって、学習モデル（予測分布） $p(x|\hat{\mu}, \hat{\sigma^2})$ が得られたことになります。
ここで、 $\hat{\mu}$ と $\hat{\sigma^2}$ はデータセットが与えられると一意に決まります。（普通にデータの平均・分散を出すだけです。）

2. 学習したモデルを基に、データの異常度合いをスコアリング（異常度の定義）

　一般に、異常度の定義として、負の対数尤度を採用されることが多いです。（これは情報理論におけるシャノン情報量と呼ばれるものです）
新たな観測値 $x'$ に対する異常度 $a(x')$ は、

$\displaystyle \begin{eqnarray} a(x') & = & -\ln p(x'|\hat{\mu}, \hat{\sigma^2}) \\ & = & -\ln \biggl[\frac{1}{\sqrt{2\pi \hat{\sigma}^2}}\exp{\left\{-\frac{(x'-\hat{\mu})^2}{2\hat{\sigma}^2}\right\}} \biggr] \\ & = & \frac{1}{2}\ln (2\pi \hat{\sigma}^2)+\frac{1}{2\hat{\sigma}^2}(x'-\hat{\mu})^2 \end{eqnarray}$

第1項は観測データ $x'$ に依存しないので無視します。全体に2をかけると、異常度は、

$\displaystyle a(x')=\biggl(\frac{x'-\hat{\mu}}{\hat{\sigma}}\biggl)^2$

という形で表されます。
分子 $(x'-\hat{\mu})$ は観測データの標本平均からのずれの大きさを表しており、ずれが大きくなると異常度が大きくなります。
分母は $\hat{\sigma}$ であり、標準偏差で割ることで、もともとばらつきの大きいものは多少のずれでも多めに見るようになっており、異常度の表現として非常に直感に近い形になっていることがわかります。

3. 閾値の設定

　異常度が決まると、それに閾値を設定することで異常判定をすることができます。閾値は経験に基づき主観的に決める方法もありますが、できるだけ客観的基準に基づいて決めるのが望ましいです。
　ホテリング理論は、観測データが正規分布に従うことを仮定した下で異常検知を行う非常に有名な古典理論であり、広く応用されています。ホテリング理論が有効な理由の一つは、異常度が従う確率分布を、明示的に導くことができる点です。異常度の確率分布が既知であれば、それに基づいて閾値を設定するのは容易です。例えば、「正常には1パーセント未満でしか起こらないくらい稀な値だから、きっと正常ではないだろう」という論理になります。つまり、パーセント値により異常判定を行うことができるようになります。
　ホテリング理論によると異常度 $a(x’)$ はデータ数Nが十分に大きい時、自由度1のカイ二乗分布に従うということが数学的に証明できます。
　下のグラフは自由度1のカイ二乗分布の確率密度関数です。

f:id:hirotsuru314:20170919224519p:plain:w550

　ここで、カイ二乗分布とは、 $Z_1,Z_2,...Z_n$ を独立な標準正規分布 $N(0, 1)$ に従う確率変数としたとき、

$\displaystyle \chi^2=Z_1^2+Z_2^2+...+Z_n^2$

とすると、確率変数 $\chi^2$ が従う確率分布を自由度 $n$ のカイ二乗分布といいます。
（異常度がカイ二乗分布に従うことの数学的な証明は難解なため避けますが、Nが十分に大きいとき、標本平均・標本分散は母平均・母分散に近づき、 $a(x’)$ が「確率変数から母平均を引いて、母標準偏差で割ったもの」すなわち標準正規分布 $N(0, 1)$ の2乗に従うようになるからだと考えてよいのだと思います）

Pythonによる実践

　これまで説明した1次元正規分布に基づく異常検知をPythonで実装したいと思います。
　ここで、入力として、データと異常度の閾値(%)を与えると、出力として、異常値とそのインデックス番号を返す関数を作りたいと思います。

import numpy as np
from scipy import stat

def hotelling_1d(data, threshold):
    """
    Parameters
    ----------
    data : Numpy array
    threshold : float

    Returns
    -------
    List of tuples where each tuple contains index number and anomalous value.
    """
    #Covert raw data into the degree of abnormality
    avg = np.average(data)
    var = np.var(data)
    data_abn = [(x - avg)**2 / var for x in data]

    #Set the threshold of abnormality
    abn_th = stats.chi2.interval(1-threshold, 1)[1]

    #Abnormality determination
    result = []
    for (index, x) in enumerate(data_abn):
        if x > abn_th:
            result.append((index, data[index]))
    return result

　この関数に引数のdataとして、200人の体重データを渡すと以下のように、

hotelling_1d(data, 0.01)
#-> [(11, 166), (20, 119)]

　異常値として2点返ってきます。
　結果を可視化すると下のようになります。縦軸は異常度です。グラフ内の赤の破線は異常度の閾値です。（プログラム中の「abn_th」の値に対応）

f:id:hirotsuru314:20170919224604p:plain:w550

　プログラムを見てもわかりますが、結局ほぼデータの平均と分散を計算して、ちょろっと異常度に変換するとできてしまいます。カイ二乗分布の確率密度関数は非常に複雑ですが、SciPyを用いるとその積分値は簡単に出せます。

ホテリング理論の課題

　ホテリング理論では、観測データ同士が互いに独立で、それぞれが単一の正規分布に従っていると仮定していますが、それは非常に強い制約と言えます。したがって、下記のようなデータには適用は困難です。

１. 複数のモードがあるような複雑な系への適用
　例えば、複数のピークを持ち、単一の正規分布で仮定できないような場合です。

２. 値が動的に変化する時系列データへの適用が困難
　正規分布を仮定するということは、データの分布がひと山でだいたい安定しているとみなしいますが、平均の値が揺らぐようなデータもあります。いわゆる時系列データと呼ばれる類のものですが、そのようなデータにホテリング理論は適用が困難です。
　時系列データの異常検知についてはまた別の機会に書きたいと思います。