ロボットと趣味と自堕落と

ロボットの事なんか一言も書いてない

機械学習による異常検知① ー正規分布に従うデータからの異常検知

C++ 技術 プログラミング 異常検知

はじめに

2年ぐらい前に下記の本を利用して異常検知をちょこっとやっていたのだけれど、気が向いたのでここに纏めていこうとおもう。 僕は確率統計苦手マンで、説明下手くそマンなので、間違っているところや説明が雑な点があったら指摘してください。*1

使う本

井手剛の『入門 機械学習による異常検知 Rによる実践ガイド』を参考に用いる。 なお、個人的にﾌﾟﾖｸﾞﾔﾐﾝｸﾞの練習でもやろうかなとか思っているので、コードはC++で書く。*2

入門 機械学習による異常検知―Rによる実践ガイド

入門 機械学習による異常検知―Rによる実践ガイド

異常検知とは

何らかの分布に従って生成されているデータから得られた規範を元に、新たなデータがその規範に従っているかどうかを判定し、規範から外れているものを異常とする仕組み。 書籍の言葉を使うと、

正常となるモデルをデータから作り、そのモデルから外れるものを異常とすること。

です。 詳しい話は本の1章に任せます。

異常検知の手順

1. 分布推定
正常なデータからモデルとなる確率分布を作り、分布のパラメータ（平均や分散など）を推定する。
2. 異常度の定義
正常からのズレ具合を定義する。例えば負の対数尤度など。
3. 閾値の設定
ホテリング理論などに基づいて、異常度に対する閾値を設定する。
4. 異常検知
定義した異常度を計算し、閾値よりも高い値を異常と判定する。

1変数正規分布に基づく異常検知

データの分布を推定する

異常な値が含まれていないとされるデータ*3 D=\left\{x_1, x_2, ..., x_N\right\} が
正規分布に従うとする。
\displaystyle \mathcal {N}( x | \mu , \sigma^2) = {\frac{1}{(2\pi\sigma^2)^{1/2}}} exp \left\{ -{\frac{1}{2\sigma^2}}(x-\mu)^2 \right\} \tag{1}
この分布のパラメータ \mu, \sigmaを最尤推定で求めると
\displaystyle \hat{\mu} = \frac{1}{N}\sum_{n=1}^{N} x_n \tag{2}

\displaystyle \hat{\sigma} = \frac{1}{N} \sum_{n=1}^{N}(x_n-\hat{\mu})^2 \tag{3}
となる。

異常度を定義

負の対数尤度を異常度 a(x)と定義する。
(1)式の対数尤度を計算すると、以下のようになる。
\displaystyle \ln{\mathcal{N}}=\frac{1}{2\hat{\sigma}^2}(x-\hat{\mu})^2+\frac{1}{2}\ln{2\pi\hat{\sigma}^2} \tag{4}
ここで xを観測値 x=x'とすると、ある観測値に対する異常度を計算する事ができる。
(4)式で観測値 x'に依存する項のみを抜き出し、2を掛けると、異常度 a(x')は
\displaystyle a(x')=\frac{1}{\hat{\sigma}^2}(x'-\hat{\mu})^2=(\frac{x'-\hat{\mu}}{\hat{\sigma}})^2 \tag{5}
と定義できる。

閾値を設定

  • ホテリング理論

1次元の観測データ D=\left\{x_1, x_2, ..., x_N\right\} の各観測値が独立に同じ分布 \mathcal{N}(\mu, \sigma^2)に従い、新たな観測値 x'も同じ分布に独立に従うとする。このとき、式(5)の a(x')の定数倍は自由度 (1, N-1)の F分布に従う。すなわち
\displaystyle \frac{N-1}{N+1}a(x')～\mathcal{F}(1, N-1) \tag{6}
特に、 N \gg 1のときは、 a(x')そのものが自由度1、スケール因子1のカイ二乗分布に従う。
\displaystyle a(x')～\chi^2(1,1) \tag{7}

よって、閾値はカイ二乗分布を元に決めれば良い。
つまり、カイ二乗分布において、 a_{th} \le a(x)の面積が \alpha\%になるときの異常度を閾値をすれば良いということ。
具体的には、例えば「閾値を \alpha=15.7\%となるように選ぶ」とすると、下図の青い領域の面積が 0.157となるように異常度の限界値を決めるという事となり、その時の異常度の閾値は a_{th}=2.0となる。 f:id:i-hako:20180301014137p:plain
このことを数式に表すと以下の積分方程式を解いて閾値 a_{th}を求める事となる。*4
$$ \displaystyle \begin{align} \alpha &= \int_{a_{th}}^{\infty} \chi^2(u|1,1) du \\ &=1-\int_0^{a_{th}} \chi^2(u|1,1) du \tag{8} \end{align} $$

実際のコード

#include <iostream>
#include <boost/math/distributions/chi_squared.hpp>
using namespace std;
using namespace boost::math;

double Uniform(){
  return ((double)rand()+1.0)/((double)RAND_MAX+2.0);
}

vector<double> GenerateData(int data_size){
  vector<double> data;
  // srand((unsigned)time(NULL));
  srand(100);
  double mu=0.0, sig=10.0;
  for(int i = 0; i < data_size; i++){
    double z = sqrt(-2.0*log(Uniform()))*sin(2.0*M_PI*Uniform());
    double val = 170.0+mu+sig*z;
    data.push_back(val);
  }
  return data;
}

double ChisqrSetCriticalValue(double dof, double los){
  chi_squared dist(dof);
  double p_val = quantile(complement(dist, los));
  return p_val;
}

double CalculateAnomaly(double mu, double sig, double val){
  return pow((val-mu), 2)/sig;
}

int main(void){
  // 正規分布に従う適当な正常データを生成
  vector<double> normal_data = GenerateData(100);

  // 分布のパラメータを推定
  double mu = 0.0, sig = 0.0;
  for(int i = 0; i < normal_data.size(); i++)
    mu += normal_data[i]/normal_data.size();
  for(int i = 0; i < normal_data.size(); i++)
    sig += pow(normal_data[i]-mu,2)/normal_data.size();

  // 異常なデータを生成して元データにくっつける
  vector<double> abnormal_data = GenerateData(3);
  for(int i = 0; i < abnormal_data.size(); i++){
    double data_value = abnormal_data[i];
    normal_data.push_back(data_value + 50);
  }

  // 閾値を計算 (α=0.01)
  double a_th = ChisqrSetCriticalValue(1, 0.01);
  cout << "a_th: " << a_th << endl;

  // 異常度を計算
  vector<double> anomaly_data;
  for(int i = 0; i < normal_data.size(); i++){
    double data_val = normal_data[i];
    double anomaly = CalculateAnomaly(mu, sig, data_val);
    if(a_th <= anomaly){
      cout << "index: " << i << " is abnormal data" << endl;
      cout << "value: " << data_val << endl;
      cout << "anomaly score: " << anomaly << endl;
    }
  }

  // 描画はお好きなように…
  // (csvで吐き出してgnuplotでもいいし、matplotでもいいと思います)
  return 0;
}

閾値を求める際に、boost::mathのquantile()を使いました。 簡単に棄却限界値を求められるので楽です。

実行結果

f:id:i-hako:20180301024757p:plain

a_th: 6.6349
index: 100 is abnormal data
value: 234.823
anomaly score: 38.1338
index: 101 is abnormal data
value: 221.681
anomaly score: 24.1974
index: 102 is abnormal data
value: 224.632
anomaly score: 27.052

適当にデータ生成してしまったので、怪しい結果のやつもいますが、最後に付け加えた3つの異常データがすべて閾値の点線よりも上に存在している事がわかります。

おわりに

すごくざっくり説明したので、これを読んで「？？？」となった人がいるかもしれません。*5
実際のホテリング理論は結構複雑な証明を2、3個挟む必要があるのですが、力尽きたので割愛しました。*6
また暇な時に書きたいと思います。
今度は多次元正規分布に基づく異常検知について書いていきたいと思います。

*1:可能な範囲で対応します

*2:雑魚ードなのはご容赦…

*3:異常値が含まれている場合は事前にデータクレンジングなどを行う必要がある

*4:実際には直接計算する必要は殆どなく、計算ライブラリなどに頼れば良い

*5:ごめんなさい…

*6:はてブロなんでこんなに数式書くの面倒なの…