pythonで「入門 機械学習による異常検知」を読む1 正規分布に従うデータからの異常検知

「入門 機械学習による異常検知」はRで書かれた異常検知の本です。
今回は実装部分をpythonで書き直していこうと思います。

個人的な勉強のメモ書きとなります。
導出や詳細などは「入門 機械学習による異常検知」を読んでいただければと思います。


### 1. 異常検知手順の流れ

0. 準備
$M$次元の観測値が$N$個手元にあると仮定する。  
データをまとめて$D$という記号で表し、この中には異常な観測値が含まれていないか、含まれていたとしてもその影響は無視できると仮定する。  
$$
D=\{\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)},\cdots,\boldsymbol{x}^{(N)}\}
$$

1. ステップ1(分布推定)
データの性質に応じた適切な確率分布のモデルを仮定する。 
パラメータを$\boldsymbol{\theta}$という記号で表す。
分布推定の問題とは、$p(\boldsymbol{x}|\boldsymbol{\theta})$における未知パラメータ$\boldsymbol{\theta}$ を、$D$から決める問題である。

2. ステップ2(異常度の定義)
観測値$\boldsymbol{x}$に対する予測分布を、$p(\boldsymbol{x}|D)$と表す。
新たな観測値$\boldsymbol{x'}$に対する異常度$a(\boldsymbol{x'})$を、予測分布に対する負の対数尤度、
$$
a(\boldsymbol{x'})=-\ln p(\boldsymbol{x}|D)
$$
で測ることを基本とする。

3. ステップ3(閾値の設定)
異常度が決まると、それに異常判定のための閾値を付すことで異常検知ができる。
ただし、異常度の確率分布を求めるのは簡単なことではなく、正常データ$D$における割合(分位点)を使うのが一般的である。

### 2. 1変数正規分布に基づく異常検知
Davisというデータを使用する。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/0b88bc2c-0a0e-e861-e4c6-07213f8539c3.png)

```python
plt.hist(davis['weight'], bins=20);
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/eaaf4d6e-cfe0-2255-da95-14ab63f52f37.png)

60付近をピークとした山となっているが、120や160あたりに異常値と思われるデータが含まれている。

ここで、正規分布と異常度について簡単に説明する。
確率変数を$x$としたとき、平均$\mu$、分散$\sigma^2$をもつ正規分布は$N(x|\mu,\sigma^2)$は、

```math
N(x|\mu,\sigma^2)=\frac{1}{(2\pi\sigma^2)^{1/2}}\exp\biggl\{-\frac{1}{2\sigma^2}(x-\mu)^2 \biggr\}
```

この平均$\mu$、分散$\sigma^2$は、データから決めるべきパラメータであり、最尤推定で決定するのが一般的である。
  
```math
\hat{\mu}=\frac{1}{N}\sum_{n=1}^Nx^{(n)}\\
\hat{\sigma}^2=\frac{1}{N}\sum_{n=1}^N(x^{(n)}-\hat{\mu})^2
```

それぞれ、**標本平均**と**標本分散**である。

負の対数尤度を異常度と考えることとして、**異常度**$a(x')$を次のように定義する。

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

異常度の確率分布の導き方として、**ホテリング統計量**が知られている。

>**定理1　(ホテリング統計量の分布(1変数))**
1次元の観測データ$D=\{x^{(1)},x^{(2)},\cdots,x^{(N)}\}$の各観測値が独立に同じ分布$N(\mu,\sigma^2)$に従い、新たな観測値$x'$も同じ分布に独立に従うとする。このとき、$a(x')$の定数倍は、自由度$(1,N-1)$のF分布に従う。
$$
\frac{N-1}{N+1}a(x')\sim F(1,N-1)
$$
特に、$N\gg 1$のときは、$a(x')$そのものが自由度1、スケール因子1のカイ二乗分布に従う。
$$
a(x')\sim \chi^2(1,1)
$$
$\frac{N-1}{N+1}a(x')$は、ホテリングの統計量(ホテリングの$T^2$)と呼ぶ。

次に手順の確認をする。
今回はカイ二乗分布の上側1%点を閾値として異常判定を行う。

>手順1
>**1. 準備**　異常が含まれていないか、含まれていたとしてもごく少数と思われるデータセットを用意する。異常判定の閾値を確率値$\alpha$で与え、カイ二乗分布の表から、異常度の閾値$a_{th}$を求めておく。
>**2. ステップ1(分布推定)**　標本平均および標本分散
$$
\hat{\mu}=\frac{1}{N}\sum_{n=1}^Nx^{(n)},\ \hat{\sigma}^2\frac{1}{N}\sum_{n=1}^N(x^{(n)}-\hat{\mu})^2
$$
を計算する。
>**3. ステップ2(異常度の判定)**　新たな観測値$x'$が得られるたび、異常度
$$
a(x')=\biggl(\frac{x'-\hat{\mu}}{\hat{\sigma}} \biggr)^2
$$
の値を計算する。
>**4. ステップ3(閾値判定)**　異常度が閾値$a_{th}$を超えたら異常と判定する。

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import chi2

# ステップ1
# カイ二乗分布による1%水準の閾値
ath = chi2.ppf(0.99, df=1)

# ステップ2
x = davis['weight']

mu = np.mean(x)
s2 = np.mean((x - mu_hat)**2)

# ステップ3
a = (x - mu)**2 / s2

# ステップ4
plt.scatter(np.arange(len(x)), a, color = ['red' if a > ath else 'blue' for a in a]);
plt.hlines(ath, xmin=0, xmax=200, linestyles='--', color='red');
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/6e1efdde-26c4-5e95-7422-523aa9a43e11.png)

### 3. 多変量正規分布に基づく異常検知
独立に同じ分布に従う$M$次元の$N$個の観測値からなるデータ$D$を考える。  
多次元正規分布は、

```math
N(\boldsymbol{x}|\boldsymbol{\mu},\Sigma)=\frac{|\Sigma|^{-1/2}}{(2\pi)^{M/2}}\exp\biggl\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\Sigma(\boldsymbol{x}-\boldsymbol{\mu}) \biggr\}
```

と書くことができる。
また、$\boldsymbol{\mu}$、$\Sigma$の最尤推定値を求めると次のようになる。

```math
\hat{\boldsymbol{\mu}}=\frac{1}{N}\sum_{n=1}^N\boldsymbol{x}^{(n)}\\
\hat{\Sigma}=\frac{1}{N}\sum_{n=1}^N(\boldsymbol{x}^{(n)}-\hat{\boldsymbol{\mu}})(\boldsymbol{x}^{(n)}-\hat{\boldsymbol{\mu}})^T
```

異常度は1変数の場合と同様に、負の対数尤度$-\ln N(\boldsymbol{x}|\boldsymbol{\mu},\Sigma)$の2倍をもとに次のように定義する。
$$
a(\boldsymbol{x'})=(\boldsymbol{x'}-\hat{\boldsymbol{\mu}})^T\hat{\Sigma}^{-1}(\boldsymbol{x'}-\hat{\boldsymbol{\mu}})
$$
これを、**マハラノビス距離**(の2乗)と呼ぶこともある。

ここで、1変数で定義したホテリング統計量を多変量に拡張する。

>**定理2　(多変数のホテリングの$T^2$理論)**  
$M$次元正規分布$N(\boldsymbol{\mu},\Sigma)$からの$N$個の独立標本$\{\boldsymbol{x}^{(1)},\cdots,\boldsymbol{x}^{(N)} \}$に基づき、標本平均$\boldsymbol{\mu}$、標本分散$\Sigma$を、
>
```math
\hat{\boldsymbol{\mu}}=\frac{1}{N}\sum_{n=1}^N\boldsymbol{x}^{(n)}\\
\hat{\Sigma}=\frac{1}{N}\sum_{n=1}^N(\boldsymbol{x}^{(n)}-\hat{\boldsymbol{\mu}})(\boldsymbol{x}^{(n)}-\hat{\boldsymbol{\mu}})^T
```
>と定義する。  
$N(\boldsymbol{\mu},\Sigma)$からの独立標本$\boldsymbol{x'}$を新たに観測したとき、以下が成立する。
>1. $\boldsymbol{x'}-\hat{\boldsymbol{\mu}}$は、平均$\boldsymbol{0}$、共分散$\frac{N+1}{N}\Sigma$の$M$次元正規分布に従う。
>2. $N\hat{\Sigma}$は、$\boldsymbol{x'}-\hat{\boldsymbol{\mu}}$と統計的に独立に、自由度$N-1$、スケール行列$\Sigma$の$M$次元ウィシャート分布に従う。
>3. $\frac{N-M}{(N+1)M}a(\boldsymbol{x'})$は、自由度$(M,N-M)$のF分布に従う。
>4. $N\gg M$の場合は、$a(\boldsymbol{x'})$は、近似的に、自由度$M$、スケール因子1のカイ二乗分布に従う。
>
```math
T^2=\frac{N-M}{(N+1)M}(\boldsymbol{x'}-\hat{\boldsymbol{\mu}})^T\hat{\Sigma}^{-1}(\boldsymbol{x'}-\hat{\boldsymbol{\mu}})
```
>は、ホテリング統計量(ホテリングの$T^2$)と呼ぶ。

変数の数$M$がよほど多くなければ$N\gg M$が成り立つことがほとんどである。
この場合、異常度$a$が、データの物理的単位や数値によらず、普遍的に、自由度$M$、スケール因子1のカイ二乗分布に従う。

```python
x = davis[['weight', 'height']]

plt.plot(x.weight, x.height, 'o');
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/284fc833-216c-8f70-dc2d-caec44d560c0.png)

```python
# 閾値の決定
ath = chi2.ppf(0.99, df=1)

# 最尤推定値の計算
mx = np.mean(x, axis=0) # 平均
Xc = x - mx
Sx = Xc.T @ Xc / len(x) # 共分散

# 異常度の計算
a = np.sum(Xc * np.matrix(np.linalg.inv(Sx)@Xc.T).T, axis=1)

plt.scatter(np.arange(len(x)), a, color = ['red' if a > ath else 'blue' for a in a]);
plt.hlines(ath, xmin=0, xmax=200, linestyles='--', color='red');
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/c9a313ec-774d-2252-33c0-39a9864a26e6.png)

### 4. マハラノビス=タグチ法
ホテリング理論で計算されるのは全系の総合的な異常度であり、個別の異常度ではなかった。
この課題に対応するものとして**マハラノビス=タグチ法(MT法)**がある。

>**手順2 (マハラノビス=タグチ法)**
>**1. 分布推定**　$D$を基に、標本平均と標本共分散を求める
>**2. 異常度の計算**　$D$の中の各標本に対し、1変数当たりのマハラノビス距離を計算する
>**3. 異常判定1**　$D$の標本が正常範囲に入るように1変数当たりのマハラノビス距離の閾値を求める
>**4. 異常判定2**　$D'$の各標本に対して、$M$変数の中からいくつかの変数を選び、その変数集合の1変数当たりの異常度を計算する

変数集合$q$に対するSN比を,

```math
SN_q=-10\log_{10}\biggl\{\frac{1}{N'}\sum_{n=1}^{N'}\frac{1}{a_q(\boldsymbol{x'}^{(n)})/M_q} \biggr\}
```

として導入する。  
ここで、$q$は変数の取捨選択パターンを区別する添字で、$M_q$はパターン$q$における変数の数、$a_q$はパターン$q$に対応して、$M_q×M_q$次元の共分散行列を使ったときの異常度である。  

今回はroadデータを用いて実装を行う。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/a0b68524-561d-7af2-4c67-d48cde866d09.png)

```python
# 変数はdrivers(運転者数)で割った値の対数をとったものを使用
x = road.apply(lambda x: x / road['drivers'])
del x['drivers']
x = np.log1p(x)

# 標本平均と標本共分散を求める
mx = np.mean(x)
Xc = x - mx
Sx = Xc.T@Xc/len(x)

# 異常値の計算(マハラノビス距離)
a = np.sum(Xc * np.matrix(np.linalg.inv(Sx)@Xc.T).T/x.shape[1], axis=1)

plt.scatter(np.arange(len(x)), a, color = ['red' if a > 1 else 'blue' for a in a]);
plt.hlines(1, xmin=0, xmax=27, linestyles='--', color='red');
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/20415c60-7147-d738-1542-92a9fdbc4322.png)

SN比は、$N'=1,M_q=1$として各変数のSN比を計算する。  
この場合のSN比は非常に簡単になり、

```math
SN_q=10\log_{10}\frac{a_q(\boldsymbol{x'})}{M_q}=10\log_{10}\frac{(x'_q-\hat{\mu}_q)^2}{\hat{\sigma}_q^2}
```

$\hat{\mu}_q$および$\hat{\sigma}_q^2$はそれぞれ、第$q$番目の変数の標本平均と標本分散である。
異常と判定されたもののうち、「Calif」についてSN比解析を行う。

```python
# 対象を選択
xc_prime = Xc.loc['Calif',:]

# SN比の計算
SN1 = 10*np.log(xc_prime**2 / np.diag(Sx))

plt.bar(x=Xc.columns, height=SN1);
plt.hlines(0, xmin=-0.5, xmax=4.5, linestyles='-', color='black');
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2358656/12ee786c-8edf-7c39-ae59-38a74cc27b21.png)


負のSN比は平均からの偏差が標準偏差よりも小さい場合に生じる。
Califの大きな異常度はほとんどすべてfuelに帰せられることがわかる。


以上となります。
### 次回
非正規データからの異常検知

### 参考文献

https://www.amazon.co.jp/%E5%85%A5%E9%96%80-%E6%A9%9F%E6%A2%B0%E5%AD%A6%E7%BF%92%E3%81%AB%E3%82%88%E3%82%8B%E7%95%B0%E5%B8%B8%E6%A4%9C%E7%9F%A5%E2%80%95R%E3%81%AB%E3%82%88%E3%82%8B%E5%AE%9F%E8%B7%B5%E3%82%AC%E3%82%A4%E3%83%89-%E4%BA%95%E6%89%8B-%E5%89%9B/dp/4339024910