yProcessingClub

すみません、許してください

多次元ホテリングのT^2法

今日も異常検知をやっていく.

今回は多次元のホテリングのT^2法だ.
前回はその準備編として多次元正規分布との和解を成し遂げたので和解していない人は是非ご確認願いたい.



以下のようにしてデータの異常検知を行っていく.
まず,データがM次元正規分布に従うと仮定する.

M次元正規分布は以下で定義される.


\displaystyle
{\mathcal N}(\boldsymbol{\mu } , \boldsymbol{\Sigma}) = \frac{1}{\sqrt{(2 \pi)^{M}|\boldsymbol{\Sigma}|}}   \exp{ \Biggl( -\frac{1}{2} ( \boldsymbol{x}-\boldsymbol{\mu} )^T  \boldsymbol{\Sigma}^{-1} ( \boldsymbol{x}-\boldsymbol{\mu} ) \Biggr) }


その異常度はデータ分布の自然対数を取って,さらに不要な係数などを消し飛ばすことで定義される.


\displaystyle
\begin{eqnarray}
\ln {\mathcal N}(\boldsymbol{\mu } , \boldsymbol{\Sigma})&=& \ln \frac{1}{\sqrt{(2 \pi)^{M}|\boldsymbol{\Sigma}|}} -\frac{1}{2} ( \boldsymbol{x}-\boldsymbol{\mu} )^T  \boldsymbol{\Sigma}^{-1} ( \boldsymbol{x}-\boldsymbol{\mu} ) 
\end{eqnarray}

右辺の第1項は\boldsymbol{x}に依存しないので取っちゃおう.第2項の係数も要らないので捨ててしまおう.
すると以下のようになる.


\displaystyle
\begin{eqnarray}
a(\boldsymbol{x})&=&  ( \boldsymbol{x}-\boldsymbol{\mu} )^T  \boldsymbol{\Sigma}^{-1} ( \boldsymbol{x}-\boldsymbol{\mu} ) 
\end{eqnarray}

異常度a(\boldsymbol{x})が無事に定義された.これはマハラノビス距離(の2乗)とも呼ばれる.


なお,上式は1次元正規分布の異常度を素直に多次元に拡張した形であることは数式を見比べれば分かるだろう.


\displaystyle
a ( x ) = \frac{(x-{\mu})^2}{{\sigma} ^ 2}


さて,なぜ異常度をこのように定義するかについて,お気持ちを述べておく.
データが正規分布に従うと仮定した.
正規分布の定義は上に示した通りで,その数式には\expが入っている.

あるデータ\boldsymbol{x}^{'}が異常であると言えるのはどのような時か.
それは\boldsymbol{x}^{'}が分布から外れている時であり,{\mathcal N}(\boldsymbol{x}^{'}|\boldsymbol{\mu } , \boldsymbol{\Sigma})が小さい場合であり,\expの中身-1/2 ( \boldsymbol{x}-\boldsymbol{\mu} )^T  \boldsymbol{\Sigma}^{-1} ( \boldsymbol{x}^{'}-\boldsymbol{\mu} ) が小さい場合であり,つまり( \boldsymbol{x}-\boldsymbol{\mu} )^T  \boldsymbol{\Sigma}^{-1} ( \boldsymbol{x}^{'}-\boldsymbol{\mu} ) が大きい場合である.
このようにして\expの中身を取り出すために自然対数を取っている(\ln \exp(x) = x)という話である.



さて,異常度a(\boldsymbol{x})の分布はどうなるだろう.
結論としては,その分布は自由度M,スケール因子1カイ二乗分布に近似的に従うのである.
こうなることは頭の良い人が考えたので結果だけ使わせてもらおう.

お気持ちとしては,異常度が小さいデータほど頻度が多く,異常度が大きくなるにつれてその頻度は減少していく感じだ.
自由度1カイ二乗分布は以下であり,直観に合っている.

f:id:Yuri-Processing-Club:20220122204757p:plain

次元が大きくなるにつれて自由度も大きくなるという感じだ.

さて,異常度の分布が分かったので,あとは閾値を決めるだけだ.
ある異常度a_{\rm th}が起きる確率\alphaを手動で決定し,a_{\rm th}を異常度の分布の確率密度関数から求めるだけだ.

f:id:Yuri-Processing-Club:20220122210610p:plain



\displaystyle
\alpha = \int_{a_{\rm th}}^{\infty} \chi^2(a | M, 1)da

確率密度関数は全区間積分すると1になる(\int_{-\infty}^{\infty} \chi^2(a | M, 1)da = 1)ので,定義域は0 \leqq a \leqq \inftyであることも考慮すると,上式は


\displaystyle
1-\alpha = \int_{0}^{a_{\rm th}} \chi^2(a | M, 1)da
とも書ける.

こう書き直したところで計算するのはシーキビであるが,例えばRではカイ二乗分布を求める関数が備わっているので,以下のように簡単に求められる.

# 1%水準の場合
th <- qchisq(0.99, M)

こうして求めたa_{\rm th}を使用して,新たなデータ\boldsymbol{x}^{'}の異常度a(\boldsymbol{x}^{'})閾値を超えたら警報を発すれば良い.


\displaystyle
\begin{eqnarray}
a(\boldsymbol{x}^{'})&=&  ( \boldsymbol{x}^{'}-\boldsymbol{\mu} )^T  \boldsymbol{\Sigma}^{-1} ( \boldsymbol{x}^{'}-\boldsymbol{\mu} ) > a_{\rm th}
\end{eqnarray}
 ? 異常:正常