『異常検知と変化検知』の6-7章（外れ値検出）要点まとめ

* 読んだ本：[異常検知と変化検知](https://www.amazon.co.jp/dp/B018K6C99U)
* 本には、コレに加えて↓が書いてある
    * 定理の証明や、使う式の導出


# 目的

後から要点を参照できるようにしておく


# 内容

## 第6章：サポートベクトルデータ記述法による異常検知

### 考え方

* 本書では、ラベルなしデータの異常検知に使う
    * $M$次元の$N$個のデータ$D = \bigl( \boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, \cdots , \boldsymbol{x}^{(N)} \bigr) $で考える
    * 「この中には異常標本が含まれていないor含まれていても圧倒的少数」とする
    * 「標準標本のほぼ全てを含む球を作り、その球に入らなかったものを異常と判定する」と考える
* ラベルありデータは別で勉強した↓
    * [「サポートベクターマシン」について](https://qiita.com/kazuya_minakuchi/items/5102f2d1d865980f279a)

まずは、中心$\boldsymbol{b}$半径$R$の球を使って「全てのデータを囲みながらできるだけ半径$R$を小さくする」と考えると、以下の最適化問題になる

```math
\min_{R^2, \boldsymbol{b}} R^2 \\
subject \; to \; || \boldsymbol{x}^{(n)} - \boldsymbol{b} ||^2 \leq R^2 \\
```

* これだと異常標本も全て含む球となるため、少しは外れる標本があっても良いという条件を加えると、以下の形になる
    * $u^{(n)} ( \geq 0 )$は、標本ごとの球からの外れ距離。球に含まれる場合0になる
    * Cは外れの許容度を示す定数（正則化定数）。大きいほど外れを少なく、小さいほど外れを多く許す

```math
\min_{R^2, \boldsymbol{b}, \boldsymbol{u}} R^2 + C \sum_{n=1}^{N} u^{(n)} \\
subject \; to \; || \boldsymbol{x}^{(n)} - \boldsymbol{b} ||^2 \leq R^2 + u^{(n)} \\
```

* ここで、ラグランジュの未定乗数法を使って変形していくと以下の最適化問題を解くことになる
    * $K(\boldsymbol{a}, \boldsymbol{b})$は$\boldsymbol{a}$と$\boldsymbol{b}$の内積を表す

```math
\min_{\boldsymbol{\alpha}} \sum_{n=1}^{N} \alpha_n K(\boldsymbol{x}^{(n)}, \boldsymbol{x}^{(n)}) - \sum_{n, n'=1}^{N} \alpha_n \alpha_{n'} K(\boldsymbol{x}^{(n)}, \boldsymbol{x}^{(n')}) \\
subject \; to \; 0 \leq \alpha_n \leq C \\
```

### 異常度

異常度$a(\boldsymbol{x'})$は「球から外れた距離」として以下で定義する

```math
\begin{align}
a(\boldsymbol{x'}) &= || \boldsymbol{x'}^{(n)} - \boldsymbol{b} ||^2 - R^2 \\
&= K(\boldsymbol{x'}, \boldsymbol{x'}) - 2 K(\boldsymbol{x'}, \boldsymbol{b}) + K(\boldsymbol{b}, \boldsymbol{b}) - R^2 \\
\end{align}
```

↑の最適化問題を解く途中で

```math
\boldsymbol{b} = \sum_{n=1}^{N} \alpha_n \boldsymbol{x}^{(n)}
```

が導出されるので

```math
R^2 = K(\boldsymbol{x'}, \boldsymbol{x'}) + \sum_{n_1, n_2 = 1}^{N} \alpha_{n_1} \alpha_{n_2} K(\boldsymbol{x}^{(n_1)}, \boldsymbol{x}^{(n_2)}) - 2 \sum_{n=1}^{N} \alpha_n K(\boldsymbol{x}^{(n)}, \boldsymbol{x}^{(n')})
```

となり、異常度も

```math
a(\boldsymbol{x'}) = K(\boldsymbol{x'}, \boldsymbol{x'}) + \sum_{n_1, n_2 = 1}^{N} \alpha_{n_1} \alpha_{n_2} K(\boldsymbol{x}^{(n_1)}, \boldsymbol{x}^{(n_2)}) - 2 \sum_{n=1}^{N} \alpha_n K(\boldsymbol{x'}, \boldsymbol{x}^{(n)}) - R^2
```

### まとめ

大事なのは以下

* 球の内部にある大多数の標本は$\alpha_n = 0$になる
    * 球の表面にある標本は$0 < \alpha_n < C$、球の外にある標本は$\alpha_n = C$となり、これらをサポートベクトル（球面を支えている=支持ベクトル=サポートベクトル）と呼ぶ。
* 最適化問題も異常度も、すべて元の標本{$\boldsymbol{x^{(n)}}$}の内積$K(\boldsymbol{x}, \boldsymbol{x'})$のみを通して依存している
    * 内積$K(\boldsymbol{x}, \boldsymbol{x'})$を、RBF（radial basis function、動径基底関数）カーネル$exp( - \sigma || \boldsymbol{x} - \boldsymbol{x'} ||^2 )$を使って置き換えると、簡単に$M$次元を∞次元に移すような非線形変換ができる


## 第7章：方向データの異常検知

### 確率分布

* データとしては、「長さ1のM次元ベクトル」がN本あると考える
* ベクトルの先端点を考えると、球面上に濃淡の模様ができる
* これは正規分布では表せないため、フォンミーゼス・フィッシャー分布を使う
    * フォンミーゼス・フィッシャー分布：球面上の正規分布のようなもの
    * [フォンミーゼスフィッシャー分布 | 高校数学の美しい物語](https://manabitimes.jp/math/1383)

フォンミーゼス・フィッシャー分布の確率分布$M(\boldsymbol{x} | \boldsymbol{\mu}, \kappa)$は以下（$\boldsymbol{\mu}$は平均方向、$\kappa$は集中度）

```math
\begin{align}
M(\boldsymbol{x} | \boldsymbol{\mu}, \kappa) &= \frac{ \kappa^{M/2-1} }{ (\pi/2)^{M/2} I_{M/2-1}(\kappa) } \exp( \kappa \boldsymbol{\mu}^T \boldsymbol{x}) \\
I_{\alpha}(\kappa) &= \frac{ 2^{-\alpha} \kappa^{\alpha} }{ \sqrt{\pi} \Gamma( \alpha + (1/2) ) } \int_{0}^{ \pi } d\phi sin^{ 2 \alpha } \phi e^{ \kappa cos \phi }  \\
\end{align}
```

平均方向については、以下のように最尤推定が可能。（集中度は、解析的には求められない）

```math
\hat{\boldsymbol{\mu}} = \frac{ \boldsymbol{m} }{ \sqrt{ \boldsymbol{m}^T \boldsymbol{m} } }, \; \boldsymbol{m} \equiv \frac{1}{N} \sum_{n=1}^{N} \boldsymbol{x}^{(n)}
```

### 異常度

確率分布の負の対数尤度をとり、キレイになるよう係数と付加定数を調整すると、異常度$a(\boldsymbol{x}')$は以下のようになる

```math
a(\boldsymbol{x}') = 1 - \boldsymbol{\mu}^T \boldsymbol{x}'
```

また、「$\kappa$が十分に大きい = aが十分に小さい」とき、近似的に以下のようにカイ２乗分布で表せる

```math
1 - \boldsymbol{\mu}^T \boldsymbol{x}' \sim \chi^2 (M-1, \frac{ 1 }{ 2 \kappa })
```

### 積率法によるカイ２乗分布の当てはめ

* 上では、まだカイ２乗分布のパラメータに未知の値である集中度$\kappa$が含まれている。
* その対処方法として、逆にデータからカイ２乗分布のパラメータ（自由度、スケール因子）を求める
* 異常度$a$が、自由度$m$、スケール因子$s$のカイ二乗分布に従うとすると、$a$と$a^2$の期待値は以下のように表せる

```math
\begin{align}
<a> &= \int_{0}^{\infty} da \; a \; \chi^2 (a | m, s ) = ms \\
<a^2> &= \int_{0}^{\infty} da \; a^2 \; \chi^2 (a | m, s ) = m (m+2) s^2 \\
\end{align}
```

これを$m, s$について解くと、以下のようになる

```math
\begin{align}
\hat{ m }_{mo} &= \frac{ 2 <a>^2 }{ <a^2> - <a>^2 } \\
\hat{ s }_{mo} &= \frac{ <a^2> - <a>^2 }{ 2 <a> } \\
\end{align}
```

データからも以下のように$a$と$a^2$の期待値は求められるので、上の式に代入して$m, s$を得られる

```math
\begin{align}
<a> & \approx \frac{1}{N} \sum_{n=1}^{N} a^{(n)} \\
<a^2> & \approx \frac{1}{N} \sum_{n=1}^{N} {a^{(n)}}^2 \\
\end{align}
```

* 注意点
    * 多くの場合、推定された自由度は実際のデータの$M-1$よりも小さくなる
    * これは、見かけ上の次元$M$が大きくても、データの主要なばらつきが一部の変数に集中していることで起こる
    * このとき、$\hat{ m }_{mo}$を有効次元と呼ぶ

### アルゴリズム

* 訓練時
    1. 訓練データから平均方向$\boldsymbol{\mu}$を求める
    2. 各標本について異常度$a(\boldsymbol{x})$を求める
    3. カイ２乗分布を当てはめ、$m$と$s$を求める
* 運用時
    1. 異常度の閾値$a_{th}$を設定する
    2. 観測データ$\boldsymbol{x}'$に対して異常度$a(\boldsymbol{x}')$を求める
    3. $\boldsymbol{x}' > a_{th}$のとき、異常と判定し警報を出す