読んだ本：異常検知と変化検知
本には、コレに加えて↓が書いてある
- 定理の証明や、使う式の導出

目的

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

内容

第4章：近傍法による異常検知

ラベルなしデータ

$M$次元の$N$個のデータ$D = \bigl( \boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, \cdots , \boldsymbol{x}^{(N)} \bigr) $
任意の位置$\boldsymbol{x'}$と、それを中心とした半径$\epsilon$の領域$V_M(\epsilon, \boldsymbol{x'})$の体積を$|V_M(\epsilon, \boldsymbol{x'})|$とする
領域$V_M(\epsilon, \boldsymbol{x'})$に、データ$D$の要素が$k$個含まれるとする
このとき、確率密度$p(\boldsymbol{x'})$は以下のようになる

p(\boldsymbol{x'}) \approx \frac{k}{N |V_M(\epsilon, \boldsymbol{x'})|}

コレを使うと、異常度$a(\boldsymbol{x'})$は

\begin{align}
a(\boldsymbol{x'}) &= - \ln p(\boldsymbol{x'}) \\
&= - \ln k + M \ln \epsilon + Const \\
\end{align}

$k$を固定して考えると、「$\epsilon$が大きい = 既存のデータを$k$個含むような領域を作るのに、必要な半径が大きい = 周りにデータが少ない = 異常度が高い」（k近傍法、k最近傍法）
$\epsilon$を固定して考えると、「$k$が小さい = 半径$\epsilon$以内に、既存のデータが少ない = 周りにデータが少ない = 異常度が高い」（$\epsilon$近傍法）

ラベルありデータ

$y=1$(異常)のときと$y=0$(正常)のときの確率密度$p(y=1 | \boldsymbol{x'})$, $p(y=0 | \boldsymbol{x'})$を求める
任意の位置$\boldsymbol{x'}$で、近い順にデータ$D$の要素を$k$個選ぶ
$k$個中、$y=1$(異常)のデータが$N^1(\boldsymbol{x'})$個、$y=0$(正常)のデータが$N^0(\boldsymbol{x'})$個あるとすると

\begin{align}
p(y=1 | \boldsymbol{x'}) &= \frac{N^1(\boldsymbol{x'})}{k} \\
p(y=0 | \boldsymbol{x'}) &= \frac{N^0(\boldsymbol{x'})}{k} \\
\end{align}

ベイズの定理を使うと

\begin{align}
p(\boldsymbol{x'} | y=1) &= \frac{p(y=1 | \boldsymbol{x'}) p(\boldsymbol{x'}) }{p(y=1)} \\
p(\boldsymbol{x'} | y=0) &= \frac{p(y=0 | \boldsymbol{x'}) p(\boldsymbol{x'}) }{p(y=0)} \\
\end{align}

よって、$p(y=1)$を$\pi_1$(全標本のうち異常の割合), $p(y=0)$を$\pi_0$(全標本のうち正常の割合)とすると、異常度$a(\boldsymbol{x'})$は

\begin{align}
a(\boldsymbol{x'}) &= \ln \frac{p(\boldsymbol{x'} | y=1)}{p(\boldsymbol{x'} | y=0)} \\
&= \ln \frac{ \pi_0 N^1(\boldsymbol{x'}) }{ \pi_1 N^0(\boldsymbol{x'}) } \\
\end{align}

マージン最大化近傍法（ラベルありデータ）

ユークリッド距離の代わりに、$M$×$M$の半正定値行列$A$を使って、2標本間の距離の2場を以下のように表す

d_A^2(\boldsymbol{x'}, \boldsymbol{x''}) = (\boldsymbol{x'} - \boldsymbol{x''})^T A (\boldsymbol{x'} - \boldsymbol{x''})

この$A$の中身を変化させて座標を伸縮・回転させ、「同じラベル同士の標本の距離は近く」「異なるラベル同士の標本の距離は遠く」なるようにする
本では、半正定値計画という最適化問題として扱い、勾配法と固有値計算を組み合わせて解いている

第5章：混合分布モデルによる逐次更新型異常検知

混合分布モデル

実際の系では、複数のモード（確率分布）が混ざったものがあり、混合分布モデルはそれを直感的に表したもの。
確率分布は以下のようになる
- $K$：系の中にあるモード（要素、コンポーネント、クラスター）の数
- $\pi_k$：系の中での、第$k$クラスターの重み（混合重み）
- $\Theta$：未知パラメータの集まり。$\Theta = (\pi_1, \cdots, \pi_K, \boldsymbol{\theta}_1, \cdots, \boldsymbol{\theta}_K)$

\begin{align}
& p(\boldsymbol{x} | \Theta) = \sum_{k=1}^{K} \pi_k p_k(\boldsymbol{x} | \boldsymbol{\theta}_k) \\
& \sum_{k=1}^{K} \pi_k = 1 \\
\end{align}

データ$D=(\boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, \cdots, \boldsymbol{x}^{(N)})$があるとき、$\Theta$についての対数尤度は

\begin{align}
L(\Theta | D) &= \ln \prod_{n=1}^{N} p(\boldsymbol{x}^{(n)} | \Theta) \\
&= \sum_{n=1}^{N} \ln \sum_{k=1}^{K} \pi_k p_k(\boldsymbol{x}^{(n)} | \boldsymbol{\theta}_k) \\
\end{align}

↑だと$N$個のデータ$D=(\boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, \cdots, \boldsymbol{x}^{(N)})$を平等に扱っているが、時刻$t$までデータ$D^{(t)}=(\boldsymbol{x}^{(1)}, \boldsymbol{x}^{(2)}, \cdots, \boldsymbol{x}^{(t)})$が得られて時刻$t$のときそれぞれのデータに$w_t^{(n)}$の重み付けをすると考えると、対数尤度は以下の形になる

L(\Theta | D^{(t)}) = \sum_{n=1}^{t} w_t^{(n)} \ln \sum_{k=1}^{K} \pi_k p_k(\boldsymbol{x}^{(n)} | \boldsymbol{\theta}_k)

また、確率分布$p_k(\boldsymbol{x} | \boldsymbol{\theta}_k)$として正規分布$N(\boldsymbol{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$を選ぶと

L(\Theta | D^{(t)}) = \sum_{n=1}^{t} w_t^{(n)} \ln \sum_{k=1}^{K} \pi_k N(\boldsymbol{x}^{(n)} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)

となる。これを最大化する$\Theta = (\pi_1, \cdots, \pi_K, \boldsymbol{\mu}_1, \cdots, \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_1, \cdots, \boldsymbol{\Sigma}_K)$を求めたい

解き方

標本のクラスターへの帰属度を導入し、「イエンセンの不等式」「EM法」を使う
- 帰属度$q_k^{(n)}$：標本$\boldsymbol{x}^{(n)}$が第$k$クラスターに所属する確率。負担率ともいう
- イエンセンの不等式：凸関数を使った不等式
- EM法(expectation-maximization algorithm)：もともとのパラメータ{$\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k$}と{$q_k^{(t)}$}の２種類の未知量を、「片方を既知として片方を求める」として両方を求めていく方法

まず、帰属度$q_k^{(n)}$を導入してイエンセンの不等式を適用した以下を最大化すると考える

\begin{align}
L(\Theta | D^{(t)}) &= \sum_{n=1}^{t} w_t^{(n)} \ln \sum_{k=1}^{K} q_k^{(n)} \frac{ \pi_k N(\boldsymbol{x}^{(n)} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) }{ q_k^{(n)} }  \\
&\geq \sum_{n=1}^{t} w_t^{(n)} \sum_{k=1}^{K} q_k^{(n)} \ln \frac{ \pi_k N(\boldsymbol{x}^{(n)} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) }{ q_k^{(n)} }  \\
\end{align}

EM法によってこれを最大化するようなパラメータを求める。（個々の計算で、ラグランジュの未定乗数法を用いる。途中の計算は割愛）
パラメータ{$\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k$}を固定し、{$q_k^{(t)}$}を求める

q_k^{(n)} = \frac{ \pi_k N(\boldsymbol{x}^{(n)} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) }{\sum_{l=1}^{K} \pi_l N(\boldsymbol{x}^{(n)} | \boldsymbol{\mu}_l, \boldsymbol{\Sigma}_l) }

パラメータ{$q_k^{(t)}$}を固定し、{$\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k$}を求める

\begin{align}
\hat{\pi}_k^{(t)} &= \frac{ 1 }{ \sum_{n'=1}^{t} w_t^{(n')} } \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} \\
\hat{\boldsymbol{\mu}}_k^{(t)} &= \frac{ 1 }{ \sum_{n'=1}^{t} w_t^{(n')} q_k^{(n')} } \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} \boldsymbol{x}^{(n)} \\
\hat{\boldsymbol{\Sigma}}_k^{(t)} &= \frac{ 1 }{ \sum_{n'=1}^{t} w_t^{(n')} q_k^{(n')} } \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} \boldsymbol{x}^{(n)} {\boldsymbol{x}^{(n)}}^{T} - \hat{\boldsymbol{\mu}}_k^{(t)} {\hat{\boldsymbol{\mu}}_k^{(t)}}^{T} \\
\end{align}

スムージングの追加

初期値の与え方が不適切だと、$\pi_k$が極端に小さい（標本がほぼ属さない）クラスターが出てきて数値的に不安定になることがある
↑を回避するため、$\pi_k$に事前分布を設定して、ゲタを履かせる（データがない状態だったら、全クラスターの混合重みが一様であるとしておく）
そうしたときの、$\hat{\pi}_k^{(t)}$の算出式は以下になる

\begin{align}
\tilde{\pi}_k^{(t)} &= \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} \\
\hat{\pi}_k^{(t)} &= \frac{ \tilde{\pi}_k^{(t)} + \gamma }{ K \gamma + \sum_{l=1}^{K} \tilde{\pi}_l^{(t)} } \\
\end{align}

逐次更新型異常検知モデル

これまでの方法だと、時刻$t$でのパラメータを求めるのに$t$個のデータ全てを使って計算する必要がある
- このままでは、データが1つ追加されるごとに全てのデータで計算しなおすことになる
データが1つ追加されたら、それまでのパラメータに追加されたデータの影響だけ反映して更新できるようにしたい
- ->逐次更新型異常検知モデル
そのために、忘却率（割引率）$\beta (0<\beta<1)$を使い、$w_t^{(n)}$を以下のように表す
- 「重み$w_t^{(n)}$が、$n$が大きい（tに近い）ほど大きく、$n$が小さい（tに近い）ほど小さく」なる。すなわち、「最近のものほど影響が多く、昔のものほど影響が小さい（忘れる）」

w_t^{(n)} = \beta (1-\beta)^{t-n}

このとき、

\begin{align}
\tilde{\pi}_k^{(t)} &= \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} \\
\tilde{\boldsymbol{\mu}}_k^{(t)} &= \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} \boldsymbol{x}^{(n)} \\
\tilde{\boldsymbol{\Sigma}}_k^{(t)} &= \sum_{n=1}^{t} w_t^{(n)} q_k^{(n)} \boldsymbol{x}^{(n)} { \boldsymbol{x}^{(n)} }^{T} \\
\end{align}

とおくと、

\begin{align}
\tilde{\pi}_k^{(t)} &= (1-\beta) \tilde{\pi}_k^{(t-1)} + \beta q_k^{(t)} \\
\tilde{\boldsymbol{\mu}}_k^{(t)} &= (1-\beta) \tilde{\boldsymbol{\mu}}_k^{(t-1)} + \beta q_k^{(t)} \boldsymbol{x}^{(t)} \\
\tilde{\boldsymbol{\Sigma}}_k^{(t)} &= (1-\beta) \tilde{\boldsymbol{\Sigma}}_k^{(t-1)} + \beta q_k^{(t)} \boldsymbol{x}^{(t)} { \boldsymbol{x}^{(t)} }^{T} \\
\end{align}

\begin{align}
\hat{\pi}_k^{(t)} &= \frac{ \tilde{\pi}_k^{(t)} + \gamma }{ K \gamma + \sum_{l=1}^{K} \tilde{\pi}_l^{(t)} } \\
\hat{\boldsymbol{\mu}}_k^{(t)} &= \frac{ \tilde{\boldsymbol{\mu}}_k^{(t)} }{ \tilde{\pi}_k^{(t)} } \\
\hat{\boldsymbol{\Sigma}}_k^{(t)} &= \frac{ \tilde{\boldsymbol{\Sigma}}_k^{(t)} }{ \tilde{\pi}_k^{(t)} } - \hat{\boldsymbol{\mu}}_k^{(t)} {\hat{\boldsymbol{\mu}}_k^{(t)}}^{T} \\
\end{align}

と表すことができ、時刻$t-1$でのパラメータ{$\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k$}に時刻$t$でのデータ$\boldsymbol{x}^{(t)}$を反映させ、時刻$t$でのパラメータを得ることができる。

パラメータが得られたら、以下の式で異常度を計算し、設定した閾値$a_{th}$を超えたら異常と判定する

a(\boldsymbol{x}) = - \ln \sum_{k=1}^{K} \pi_k N(\boldsymbol{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)

個人的な疑問点

パラメータ{$\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k$}の初期値はどう設定する？（ドメイン知識から決める？）
時刻1ステップごとに収束計算は必要ない？
- 時刻$t-1$のパラメータから計算した、時刻$t$の$q_k^{(t)}$は妥当？
- $t-1$時点のパラメータが妥当なもので、1ステップでは大きくは変わらないという想定？その場合は初期値が大事な気がする
↑の解決策（想像）
- $t-1$から$t$の逐次更新方法を使う前に、$t-1$個のデータでEM法を使った収束計算を回しておく（時刻のステップは進めない）
- $t-1$時点で妥当そうなパラメータを得てから、逐次更新を使って時刻ステップを進めていく
こうすると、「適当に決めた初期値から$t-1$時点のデータに合うパラメータを求める」「$t-1$時点で妥当なパラメータがある状態で、次の時刻$t$に進む」とできそう