ホテリングT2法による異常検知

# 内容

標準正規分布とカイ二乗分布の関係を用いて、ホテリングT2法という簡単な異常検知モデルを作成することができます。

確率変数$X$が標準正規分布$N(0,1)$に従う時、$X^2$は自由度1のカイ二乗分布$\chi^2_1$に従う為、この分布をベースに「異常である」の閾値$a$を設定し、新たに観測した値から導かれる異常度$a(x')$が閾値$a$を超えるなら「異常」とする方法です。

<img width="800" alt="065c9039-84d6-47f7-8bbd-3389bb5160ae.png名" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/635079/cb19d0c6-21da-916b-1a13-abc1a446d396.png">

実際に標準正規分布から取得した1000の標本を2乗して、プロットしてみます。

```python
from scipy.stats import norm chi2
import matplotlib.pyplot as plt
import seaborn as sns

# 実際にやってみる
x = norm.rvs(loc=0, scale=1, size=1000, random_state=100)
plt.figure(figsize=(6, 6))
sns.histplot(x)
plt.title('$N(0, 1)$からの標本 (n=1000)')
plt.savefig('正規分布標本.png')
plt.show()
```

1000の標本をそのまま取得したものをプロットします。
![3628a114-adbe-4359-bfa8-c981e9cd9961.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/635079/33b570c0-8ca3-5966-a891-0e763a0b58c1.png)

2乗してプロットします。

```python
# 二乗する
x2 = x**2
plt.figure(figsize=(6, 6))
sns.histplot(x2)
plt.title('$N(0, 1)$標本を2乗')
plt.show()
```

自由度1のカイ二乗分布に従っています（証明は[こちら](#証明)）

![標本の2乗結果.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/635079/995f2343-7e1e-26df-284f-57206376d126.png)



# 実装
先に閾値$a$を決めておきます。今回は自由度1のカイ二乗分布の97.5%点である**5.02**を閾値$a$とします。

次に、正常な標本がほとんどだと言えるデータから標本平均と標本標準偏差を求めます。
今回は例として、正常な標本が$N(10, 5^2)$に従っているとして、サンプルを取得します。

```python
# 正常時のサンプルを生成する
mu = 10
scale= 5
sample = norm.rvs(loc=mu, scale=scale, size=1000, random_state=46)

# プロット
plt.figure(figsize=(6, 6))
sns.histplot(sample)
plt.title('正常時サンプルのヒストグラム')
plt.show()
```


![正常時サンプルのヒストグラム.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/635079/8c7a1831-5707-01ca-a88a-97356c3f9aed.png)


標本平均と標本標準偏差を求めます。

```python
# 標本から平均と標準偏差を算出
sample_mu = sample.mean()
sample_sigma = sample.std()
print(sample_mu, sample_sigma)
```

```
9.92654000756153 5.092873537096135
```


新たな観測値を得たら、平均と標準偏差を用いて値を標準化し、標準化した値を2乗して異常度とします。

```python
# 新しく観測したデータ
new_sample = np.array([14, 0, 30])

# 標準化した値を2乗する
standardization_value = ((new_sample - sample_mu) / sample_sigma)
outlier_check_value = standardization_value ** 2
outlier_check_value
```

```
array([ 0.63973649,  3.79900627, 15.53526583])
```

3つの異常度$a(x')$が求まりました。
閾値$a$より大きい、3番目の観測値**30**を異常とします。


手軽に実装できますが、論理的に予測される値と実際が食い違うことから「異常」とされやすいという欠点もあります。
モーメント法を用いたカイ二乗分布の当て嵌めを行い、自由度を決定する方法もあります。


# 証明
確率変数$Z$が標準正規分布$N(0,1)$に従う時、$Z^2$が自由度$1$のカイ二乗分布$\chi^2_1$に従うことを証明します。


自由度$n$のカイ二乗分布$\chi^2_n$の確率密度関数は以下の式で与えられます。

```math
f(x) = \frac{1}{\Gamma(n/2)} \Bigl( \frac{1}{2} \Bigr)^{n/2}x^{\frac{n}{2}-1}\exp(-x/2), \quad x > 0
```

よって、自由度1の場合は、以下となります。

```math
f(x) = \frac{1}{\Gamma(1/2)} \Bigl( \frac{1}{2} \Bigr)^{1/2}x^{-\frac{1}{2}}\exp(-x/2), \quad x > 0
```


## 変数変換
$N(0,1)$に従う確率変数の二乗が$\chi^2_1$になるまでの導出について、変数変換を用います。

関数$g(\cdot)$を通して確率変数$X$を$Y = g(X)$に変換したとき、$Y$の分布を$X$の分布から導くことを考えます。
$Y$の累積分布関数$F(y)$は

```math
F(y) = P(Y \leq y)
```

であり、$Y = g(X)$より

```math
F(y) = P(g(X) \leq y) = P(X \in \{x | g(x) \leq y \})
```

となり、確率密度関数は累積分布関数の微分で得られることから、

```math
f(y) = \frac{d}{dy}F(y) = \frac{d}{dy}P(X \in \{x | g(x) \leq y \})
```

と導かれます。


### 平方変換
確率変数$X$の確率密度関数を$f(x)$としたときの、$X$の平方変換$Y = X^2$について考えます。

変数変換式から、$y > 0$に対して

```math
\{ x | x^2 \leq y \} = \{ x | - \sqrt{y} \leq x \leq \sqrt{y} \}
```

であるので、確率密度関数$g(y)$は

```math
g(y) = \frac{d}{dy}G(y) = \frac{d}{dy}P(X \in \{ x |  - \sqrt{y} \leq x \leq \sqrt{y}  \}) \\
= \frac{d}{dy}\int_{-\sqrt{y}}^{\sqrt{y}}f(x) dx \\
= \{ f(\sqrt{y}) + f(-\sqrt{y}) \}\frac{1}{2 \sqrt{y}}
```


となることから

```math
g(y) = \{ f(\sqrt{y}) + f(-\sqrt{y}) \}\frac{1}{2 \sqrt{y}}
```

と導かれます。
確率密度関数$f(x)$が$y$軸に対して対称（正規分布等）な場合は$f(-x) = f(x)$となるため、

```math
g(y) = \frac{f(\sqrt{y})}{\sqrt{y}}
```

とも書けます。

## 平方変換の適用

確率変数$X$の確率密度関数$f(x)$を$N(0,1)$として変換します。
正規分布$N(\mu, \sigma^2)$の確率密度関数は以下の式で与えられます。

```math
f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \biggl(-\frac{(x - \mu)^2}{2\sigma^2} \biggr)
```

よって、標準正規分布$N(0, 1)$は以下となります。

```math
f(x) = \frac{1}{\sqrt{2 \pi}} \exp \biggl(-\frac{x^2}{2} \biggr)
```


$X$の平方変換$Y = X^2$の$g(y) = \frac{f(\sqrt{y})}{\sqrt{y}}$を適用すると、

```math
g(y) = \frac{1}{\sqrt{y}}\cdot\frac{1}{\sqrt{2 \pi}} \exp \biggl(-\frac{y}{2} \biggr) \\
= \frac{1}{\sqrt{\pi}} \Bigl( \frac{1}{2}\Bigr)^{1/2} y^{-\frac{1}{2}} \exp (-y/2), \quad y > 0
``` 

となります。


## ガンマ関数

```math
\frac{1}{\sqrt{\pi}} \Bigl( \frac{1}{2}\Bigr)^{1/2} y^{-\frac{1}{2}} \exp (-y/2)
```

上式の$\sqrt{\pi}$を$\Gamma{(1/2)}$へ変換します。

ガンマ関数$\Gamma{(x)}$は正の実数$x$に対して、以下の式で与えられます。

```math
\Gamma{(x)} = \int_{0}^{\infty} t^{x-1} e^{-t} dt
```

$\Gamma{(1/2)}$の場合、$t = u^2$の置換積分を用いて

```math
\Gamma{(1/2)} = \int_{0}^{\infty}t^{-\frac{1}{2}}e^{-t} dt \\
= \int_{0}^{\infty} u^{-1} e^{-u^2}2u du \\
= 2 \int_{0}^{\infty}e^{-u^2}du
```

となります。

### ガウス積分

$\int_{0}^{\infty}e^{-u^2}du$の部分ですが、以下のガウス積分を用いて$\pi$へ変換します。


```math
\int_{-\infty}^{\infty}e^{-u^2}du = \sqrt{\pi}
```

ガウス積分の証明は、2次元極座標$(\rho, \phi)$を用いて可能です。
ガウス積分の値を$I$と置きます。
```math
I = \int_{-\infty}^{\infty}e^{-x^2}dx \\
```

これを2乗し$I^2$とします。こちらを解いていきます。
$x = \rho \cos{\phi}, y= \rho \sin{\phi}$と置換し、$\rho = \sqrt{x^2 + y^2}$となる事も用いると、
```math
I^2 = \int_{-\infty}^{\infty}e^{-x^2}dx \int_{-\infty}^{\infty}e^{-y^2}dy \\
= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-(x^2 + y^2)}dxdy \\
= \int_{0}^{2\pi} d\phi \int_{0}^{\infty}e^{-\rho^2} \rho d\rho \\
= 2 \pi \int_{0}^{\infty}e^{-\rho^2} \rho d \rho \\
= 2 \pi \biggl[ -\frac{e^{-\rho^2}}{2}\biggr]_{0}^{\infty} \\
= \pi
```

となりますので、
```math
I = \int_{-\infty}^{\infty}e^{-x^2}dx = \sqrt{\pi}
```


よって、

```math
\Gamma{(1/2)} = 2 \int_{0}^{\infty}e^{-u^2}du = 2 \cdot \frac{\sqrt{\pi}}{2} = \sqrt{\pi}
```

となるため、$\sqrt{\pi}$は$\Gamma{(1/2)}$へ変換でき、

```math
\frac{1}{\Gamma{(1/2)}} \Bigl( \frac{1}{2}\Bigr)^{1/2} y^{-\frac{1}{2}} \exp (-y/2), \quad y > 0
```

となります。

平方変換を用いて、標準正規分布$N(0,1)$に従う確率変数の二乗が自由度11のカイ二乗分布$\chi^2_1$になる事が証明できました。


# 参考
- [現代数理統計学の基礎](https://www.kyoritsu-pub.co.jp/book/b10003681.html)
- [微分積分](https://www.iwanami.co.jp/book/b482316.html)
- [高校数学の美しい物語](https://manabitimes.jp/math)
- [異常検知と変化検知](https://www.kspub.co.jp/book/detail/1529083.html)