異常検知の基礎

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でお手軽・異常検知 [ホテリング理論編]
【統計学】正規分布とカイ二乗分布の関係を可視化してみる。

その他の記事

Other Articles

2019/09/20
Kabuku Connect 即時見積機能のバックエンド開発

2019/08/14
Maker Faire Tokyo 2019でARゲームを出展しました

2019/07/25
夏休みだョ!WebAssembly Proposal全員集合!!

2019/07/08
鵜呑みにしないで! —— 書籍『クリーンアーキテクチャ』所感 ≪null 篇≫

2019/07/03
W3C Workshop on Web Games参加レポート

2019/06/28
TypeScriptでObject.assign()に正しい型をつける

2019/06/25
カブクエンジニア開発合宿に行ってきました 2019夏

2019/06/21
Hola! KubeCon Europe 2019の参加レポート

2019/06/19
Clean Resume きれいな環境できれいな履歴書を作成する

2019/05/20
[Web フロントエンド] 状態更新ロジックをフレームワークから独立させる

2019/04/16
C++のenable_shared_from_thisを使う

2019/04/12
OpenAPI 3 ファーストな Web アプリケーション開発(Python で API 編)

2019/04/08
WebGLでレイマーチングを使ったCSGを実現する

2019/04/02
『エンジニア採用最前線』に感化されて2週間でエンジニア主導の求人票更新フローを構築した話

2019/03/29
その1 Jetson TX2でk3s(枯山水)を動かしてみた

2019/03/27
任意のブラウザ上でJestで書いたテストを実行する

2019/02/08
TypeScript で “radian” と “degree” を間違えないようにする

2019/02/05
Python3でGoogle Cloud ML Engineをローカルで動作する方法

2019/01/18
SIGGRAPH Asia 2018 参加レポート

2019/01/08
お正月だョ!ECMAScript Proposal全員集合!!

2019/01/08
カブクエンジニア開発合宿に行ってきました 2018秋

2018/12/25
OpenAPI 3 ファーストな Web アプリケーション開発(環境編)

2018/12/23
いまMLKitカスタムモデル(TF Lite)は使えるのか

2018/12/21
[IoT] Docker on JetsonでMQTTを使ってCloud IoT Coreと通信する

2018/12/11
TypeScriptで実現する型安全な多言語対応(Angularを例に)

2018/12/05
GASでCompute Engineの時間に応じた自動停止/起動ツールを作成する 〜GASで簡単に好きなGoogle APIを叩く方法〜

2018/12/02
single quotes な Black を vendoring して packaging

2018/11/14
3次元データに2次元データの深層学習の技術(Inception V3, ResNet)を適用

2018/11/04
Node Knockout 2018 に参戦しました

2018/10/24
SIGGRAPH 2018参加レポート-後編(VR/AR)

2018/10/11
Angular 4アプリケーションをAngular 6に移行する

2018/10/05
SIGGRAPH 2018参加レポート-特別編(VR@50)

2018/10/03
Three.jsでVRしたい

2018/10/02
SIGGRAPH 2018参加レポート-前編

2018/09/27
ズーム可能なSVGを実装する方法の解説

2018/09/25
Kerasを用いた複数入力モデル精度向上のためのTips

2018/09/21
競技プログラミングの勉強会を開催している話

2018/09/19
Ladder Netwoksによる半教師あり学習

2018/08/10
「Maker Faire Tokyo 2018」に出展しました

2018/08/02
Kerasを用いた複数時系列データを1つの深層学習モデルで学習させる方法

2018/07/26
Apollo GraphQLでWebサービスを開発してわかったこと

2018/07/19
【深層学習】時系列データに対する1次元畳み込み層の出力を可視化

2018/07/11
きたない requirements.txt から Pipenv への移行

2018/06/26
CSS Houdiniを味見する

2018/06/25
不確実性を考慮した時系列データ予測

2018/06/20
Google Colaboratory を自分のマシンで走らせる

2018/06/18
Go言語でWebAssembly

2018/06/15
カブクエンジニア開発合宿に行ってきました 2018春

2018/06/08
2018 年の tree shaking

2018/06/07
隠れマルコフモデル 入門

2018/05/30
DASKによる探索的データ分析(EDA)

2018/05/10
TensorFlowをソースからビルドする方法とその効果

2018/04/23
EGLとOpenGLを使用するコードのビルド方法〜libGLからlibOpenGLへ

2018/04/23
技術書典4にサークル参加してきました

2018/04/13
Python で Cura をバッチ実行するためには

2018/04/04
ARCoreで3Dプリント風エフェクトを実現する〜呪文による積層造形映像制作の舞台裏〜

2018/04/02
深層学習を用いた時系列データにおける異常検知

2018/04/01
音声ユーザーインターフェースを用いた新方式積層造形装置の提案

2018/03/31
Container builderでコンテナイメージをBuildしてSlackで結果を受け取る開発スタイルが捗る

2018/03/23
ngUpgrade を使って AngularJS から Angular に移行

2018/03/14
Three.jsのパフォーマンスTips

2018/02/14
C++17の新機能を試す〜その1「3次元版hypot」

2018/01/17
時系列データにおける異常検知

2018/01/09
three.ar.jsを使ったスマホAR入門

2017/12/17
Python OpenAPIライブラリ bravado-core の発展的な使い方

2017/12/15
WebAssembly(wat)を手書きする

2017/12/14
AngularJS を Angular に移行: ng-annotate 相当の機能を TypeScrpt ファイルに適用

2017/12/08
Android Thingsで4足ロボットを作る ~ Android ThingsとPCA9685でサーボ制御)

2017/12/06
Raspberry PIとDialogflow & Google Cloud Platformを利用した、3Dプリンターボット(仮)の開発 (概要編)

2017/11/20
カブクエンジニア開発合宿に行ってきました 2017秋

2017/10/19
Android Thingsを使って3Dプリント戦車を作ろう ① ハードウェア準備編

2017/10/13
第2回 魁!! GPUクラスタ on GKE ~PodからGPUを使う編~

2017/10/05
第1回 魁!! GPUクラスタ on GKE ~GPUクラスタ構築編~

2017/09/13
「Maker Faire Tokyo 2017」に出展しました。

2017/09/11
PyConJP2017に参加しました

2017/09/08
bravado-coreによるOpenAPIを利用したPythonアプリケーション開発

2017/08/23
OpenAPIのご紹介

2017/08/18
EuroPython2017で2名登壇しました。

2017/07/26
3DプリンターでLチカ

2017/07/03
Three.js r86で何が変わったのか

2017/06/21
3次元データへの深層学習の適用

2017/06/01
カブクエンジニア開発合宿に行ってきました 2017春

2017/05/08
Three.js r85で何が変わったのか

2017/04/10
GCPのGPUインスタンスでレンダリングを高速化

2017/02/07
Three.js r84で何が変わったのか

2017/01/27
Google App EngineのFlexible EnvironmentにTmpfsを導入する

2016/12/21
Three.js r83で何が変わったのか

2016/12/02
Three.jsでのクリッピング平面の利用

2016/11/08
Three.js r82で何が変わったのか

2016/12/17
SIGGRAPH 2016 レポート

2016/11/02
カブクエンジニア開発合宿に行ってきました 2016秋

2016/10/28
PyConJP2016 行きました

2016/10/17
EuroPython2016で登壇しました

2016/10/13
Angular 2.0.0ファイナルへのアップグレード

2016/10/04
Three.js r81で何が変わったのか

2016/09/14
カブクのエンジニアインターンシッププログラムについての詩

2016/09/05
カブクのエンジニアインターンとして3ヶ月でやった事 〜高橋知成の場合〜

2016/08/30
Three.js r80で何が変わったのか

2016/07/15
Three.js r79で何が変わったのか

2016/06/02
Vulkanを試してみた

2016/05/20
MakerGoの作り方

2016/05/08
TensorFlow on DockerでGPUを使えるようにする方法

2016/04/27
Blenderの3DデータをMinecraftに送りこむ

2016/04/20
Tensorflowを使ったDeep LearningにおけるGPU性能調査

→
←

関連職種

Recruit

フロントエンドエンジニア(TypeScript)

業務内容

当ポジションは弊社Webサービスのフロントエンド機能設計及び実装を担当します。 サービス毎の開発チームで2週間スプリントのスクラム開発を実施しています。 週次で開発チームミーティングを実施し、実装設計の相談や工数見積もりを行います。 全ての開発コードはレビューと自動テストによって品質を保っています。 また、リファクタリングやフレームワークのバージョンアップも開発フローに組込み、技術的負債を放置しない開発を目指しています。

インターン(Webエンジニア)

業務内容

業務から独立した、調査・研究系のタスクをおまかせしています。コードレビュー、 社内での報告会、 ブログ記事執筆を通して着実にスキルアップしていただくことを目指しています。 (希望があれば、プロダクトの開発業務もおまかせします。)

→
←

お客様のご要望に「Kabuku」はお応えいたします。
ぜひお気軽にご相談ください。

お電話でも受け付けております
03-6380-2750
営業時間:09:30~18:00
※土日祝は除く