異常検知の基礎
異常検知の基礎
はじめに
カブクで機械学習エンジニアをしている大串正矢です。今回は異常検知の基礎について書きます。
背景
異常検知の手法は多種多様に存在していますがウェブ上にまとまった情報が日本語でないため記述することにしました。ただ全ての内容をこの記事で記述すると長くなるため今回は基礎的な内容にフォーカスして記述します。この記事の内容をベースに他の異常検知に関する記事を記述する予定です。
異常検知の基本的なアプローチ
下記の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)
# 生成した乱数を２乗する
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
2022/06/03
拡張子に Web アプリを関連付ける File Handling API の使い方
2022/03/22<selectmenu> タグできる子; <select> に代わるカスタマイズ可能なドロップダウンリスト
2022/03/02
Java 15 のテキストブロックを横目に C# 11 の生文字列リテラルを眺めて ECMAScript String dedent プロポーザルを想う
2021/10/13
Angularによる開発をできるだけ型安全にするためのKabukuでの取り組み
2021/09/30
さようなら、Node.js
2021/09/30
Union 型を含むオブジェクト型を代入するときに遭遇しうるTypeScript型チェックの制限について
2021/09/16
[ECMAScript] Pipe operator 論争まとめ – F# か Hack か両方か
2021/07/05
TypeScript v4.3 の機能を使って immutable ライブラリの型付けを頑張る
2021/06/25
Denoでwasmを動かすだけの話
2021/05/18
DOMMatrix: 2D / 3D 変形（アフィン変換）の行列を扱う DOM API
2021/03/29
GoのWASMがライブラリではなくアプリケーションであること
2021/03/26
Pythonプロジェクトの共通のひな形を作る
2021/03/25
インラインスタイルと Tailwind CSS と Tailwind CSS 入力補助ライブラリと Tailwind CSS in JS
2021/03/23
Serverless NEGを使ってApp Engineにカスタムドメインをワイルドカードマッピング
2021/01/07
esbuild の機能が足りないならプラグインを自作すればいいじゃない
2020/08/26
TypeScriptで関数の部分型を理解しよう
2020/06/16
[Web フロントエンド] esbuild が爆速すぎて webpack / Rollup にはもう戻れない
2020/03/19
[Web フロントエンド] Elm に心折れ Mint に癒しを求める
2020/02/28
さようなら、TypeScript enum
2020/02/14
受付のLooking Glassに加えたひと工夫
2020/01/28
カブクエンジニア開発合宿に行ってきました 2020冬
2020/01/30
Renovateで依存ライブラリをリノベーションしよう 〜 Bitbucket編 〜
2019/12/27
Cloud Tasks でも deferred ライブラリが使いたい
2019/12/25
*, ::before, ::after { flex: none; }
2019/12/21
Top-level awaitとDual Package Hazard
2019/12/20
Three.jsからWebGLまで行きて帰りし物語
2019/12/18
Three.jsに入門＋手を検出してAR.jsと組み合わせてみた
2019/12/04
WebXR AR Paint その２
2019/11/06
GraphQLの入門書を翻訳しました
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/03/29
その1 Jetson TX2でk3s(枯山水)を動かしてみた
2019/04/02
『エンジニア採用最前線』に感化されて2週間でエンジニア主導の求人票更新フローを構築した話
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
【深層学習】時系列データに対する１次元畳み込み層の出力を可視化
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