内容

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

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

065c9039-84d6-47f7-8bbd-3389bb5160ae.png名

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

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の標本をそのまま取得したものをプロットします。

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

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

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

実装

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

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

# 正常時のサンプルを生成する
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()

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

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

9.92654000756153 5.092873537096135

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

# 新しく観測したデータ
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$の確率密度関数は以下の式で与えられます。

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の場合は、以下となります。

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)$は

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

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

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

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

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$に対して

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

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

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}}

となることから

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

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

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

とも書けます。

平方変換の適用

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

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

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

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}}$を適用すると、

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

となります。

ガンマ関数

\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$に対して、以下の式で与えられます。

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

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

\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$へ変換します。

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

ガウス積分の証明は、2次元極座標$(\rho, \phi)$を用いて可能です。
ガウス積分の値を$I$と置きます。

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}$となる事も用いると、

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

となりますので、

I = \int_{-\infty}^{\infty}e^{-x^2}dx = \sqrt{\pi}

よって、

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

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

\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$になる事が証明できました。

参考