混合ガウス分布の最尤推定を解く手法のひとつとしてEMアルゴリズムを導入する．
混合ガウス分布
$K$ 個のクラスがあり，それぞれを $z_k$ と表記する．そしてその実現値を1-of-K符号化で $\boldsymbol{z}$ と表す(どれか一つだけが1，残りが0)．各クラスの周辺分布は $p(\boldsymbol{z}) = \pi_k$ であるとする．そしてクラスが既知である場合の条件付き確率は
$$\begin{align*} p (\boldsymbol{x} \mid z_{k}=1) = \mathcal{N} (\boldsymbol{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}) \end{align*}$$
であるとすると， $\boldsymbol{x}$ の周辺分布は
$$\begin{align*} p(\boldsymbol{x})=\sum_{k=1}^{K} \pi_{k} \mathcal{N} (\boldsymbol{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}) \end{align*}$$
となる．逆に $\boldsymbol{x}$ が与えられた場合の $\boldsymbol{z}$ の条件付き確率は
$$\begin{align*} \gamma (z_k) = p(z_k = 1 \mid \boldsymbol{z}) &= \dfrac{p(z_k = 1)p(\boldsymbol{x} \mid z_k = 1)}{\sum_{j=1}^{K}p(z_j = 1)p(\boldsymbol{x} \mid z_j = 1)} \\ &= \dfrac{\pi_{k} \mathcal{N} (\boldsymbol{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}(\boldsymbol{x} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j})} \end{align*}$$
で与えられる．これは $\boldsymbol{x}$ がクラスkを説明する度合いであると言え，負荷率と呼ばれる．
最尤推定の特異性
GMMから得られた観測値の集合 $\{ \boldsymbol{x}_1, \ldots, \boldsymbol{x}_N \}$ について $\mathbf{X}$ を各行が $\boldsymbol{x}_{i}^{\text{T}}$ である行列を $\mathbf{X}$ とする．すると尤度関数の対数は
$$\begin{align*} \ln p(\mathbf{X} \mid \boldsymbol{\pi}, \boldsymbol{\mu}, \mathbf{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N} (\boldsymbol{x}_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k})\right\} \end{align*}$$
で与えられる．ここではシンプルに $\mathbf{\Sigma}_k = \sigma_{k}^{2} \mathbf{I}$ とし，訓練データのうちの一つ $\boldsymbol{x}_n$ がある要素の平均 $\boldsymbol{\mu}_j$ と等しかったとする．このとき指数関数部が1になるので上の式の $\ln$ の和の中に
$$\begin{align*} \mathcal{N}(\mathbf{x}_{n} \mid \mathbf{x}_{n}, \sigma_{j}^{2} \mathbf{I})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{\sigma_{j}^{D}} \end{align*}$$
の項が現れる．補足 で見るように指数がかかった項は $\sigma \rightarrow 0$ で0に収束するが，この項は無限大に発散してしまう．よってこの対数尤度関数は無限大に発散し，尤度関数も発散する．
補足
1次元のガウス分布に従う $x \sim \mathcal{N}(\mu, \sigma^2)$ から得られたデータ $[x_1, \ldots, x_N]$ のうち $x_i$ が $\mu$ に一致しているとする．これに対して尤度関数を求めると
$$\begin{align*} \frac{1}{\sqrt{2\pi} \sigma} \prod_{n \neq i}^N \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{(x_n-\mu)^2}{2\sigma^2}} \end{align*}$$
となる．ここで $\sigma \rightarrow 0$ とすると指数部の収束が $1 / \sigma$ の発散より強いため0に収束する．
GMMのEMアルゴリズム
対数尤度を $\boldsymbol{\mu}_k$ について微分することで以下を得る．
$$\begin{align*} \gamma_{nk} &= \dfrac{\pi_k \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k)}{\sum_j \pi_k \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_j, \mathbf{\Sigma}_j)} \\ N_k &= \sum_{n=1}^{N} \gamma_{nk} \end{align*}$$
さらにこの値を使って
$$\begin{align*} \overline{\boldsymbol{\mu}_k} &= \dfrac{1}{N_k} \sum_{n=1}^{N} \gamma_{nk} \boldsymbol{x}_n \\ \overline{\mathbf{\Sigma}_{k}} &= \dfrac{1}{N_{k}} \sum_{n=1}^{N} \gamma_{nk} (\boldsymbol{x}_{n} - \overline{\boldsymbol{\mu}_{k}}) (\boldsymbol{x}_{n} - \overline{\boldsymbol{\mu}_{k}})^{\text{T}} \\ \overline{\pi_k} &= \dfrac{N_k}{N} \end{align*}$$
を得る．このように得られた値で $(\boldsymbol{\mu}_k, \mathbf{\Sigma}_k, \pi_k)$ を更新する処理を繰り返す．
EMアルゴリズムの一般形
観測値 $\mathbf{X} = \{ \boldsymbol{x}_1, \ldots, \boldsymbol{x}_N \}$ が潜在変数(例えばどのGMMにおいてどのクラスに属するか，など)を含む確率分布 $p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta} )$ に従うとする．これは人間が決めるものである．そして $p( \mathbf{X} \mid \boldsymbol{\theta} )$ を最大化するような $\boldsymbol{\theta}$ を求める(ここで $\mathbf{X}$ は観測値なので定数，よって $\boldsymbol{\theta}$ の関数になる)．
EMアルゴリズムでは潜在変数の事後分布 $p( \mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta} )$ を用いる．これも以下のようにして
$$\begin{align*} p( \mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta} ) = \frac{p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{\int p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}) d \mathbf{Z}} \end{align*}$$
前提で仮定した同時分布 $p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta} )$ のみから求められる．初期推定値 $\boldsymbol{\theta}_0$ を用いて以下の処理を繰り返す．
- $p( \mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}_i )$ を計算する
- $Q(\boldsymbol{\theta})$ を $\int p (\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta} ) \ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}_i) d\mathbf{Z}$ とする
- この最大値を与える $\boldsymbol{\theta}$ を $\boldsymbol{\theta}_{i+1}$ とする
GMMのEMアルゴリズム では負荷率 $\gamma_{nk}$ を求めた部分が(1)の事後分布を求める部分に相当する．この枠組みに当てはめると $p( \mathbf{Z} \mid \boldsymbol{x}_n, \boldsymbol{\theta}_{\mathrm{old}} )$ は
$$\begin{align*} p( \mathbf{Z} \mid \boldsymbol{x}_n, \boldsymbol{\theta} ) = \prod_{k=1}^{K} \left( \dfrac{\pi_k \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_j, \mathbf{\Sigma}_j)} \right)^{z_{nk}} \end{align*}$$
で与えられる．同様に $p(\boldsymbol{x}_n, \mathbf{Z} \mid \boldsymbol{\theta})$ は
$$\begin{align*} p( \boldsymbol{x}_n, \mathbf{Z}) = \prod_{k=1}^{K} \pi_{k}^{z_{nk}}\mathcal{N}( \boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k)^{z_{nk}} \end{align*}$$
で与えられるから対数は
$$\begin{align*} \ln p( \boldsymbol{x}_n, \mathbf{Z}) = \sum_{k=1}^{K} {z_{nk}} \ln \pi_{k} + {z_{nk}} \ln \mathcal{N}( \boldsymbol{x}_n \mid \boldsymbol{\mu}, \mathbf{\Sigma}_k) \end{align*}$$
となる．ここで $\mathbf{Z}$ は1-of-K符号化により $\{ (1, 0, \ldots ,0), (0, 1, 0, \ldots, 0), \ldots, (0, 0, \ldots, 1) \}$ を動くとして
$$\begin{align*} \int p (\mathbf{Z} \mid \boldsymbol{x}_n, \boldsymbol{\theta}_{\mathrm{old}} ) \ln p(\boldsymbol{x}_n, \mathbf{Z} \mid \boldsymbol{\theta}) d\mathbf{Z} = \sum_{k=1}^{K} \gamma_{nk} \{ \ln \pi_k + \ln \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k) \} \end{align*}$$
となる．これを $\boldsymbol{x}_n$ につき和をとったのが(2)の $Q(\boldsymbol{\theta})$ に相当する． $(\boldsymbol{\mu}_k, \mathbf{\Sigma}_k, \pi_k)$ について微分した連立方程式を0にする値を求めるとGMMのEMアルゴリズム で求めた式が得られる．
そもそも潜在変数とは
個人的には説明変数として用いている $\boldsymbol{x}$ 以外の切り捨てられた変数全てのことなのかなと思った．実際説明変数を列挙しようとすれば，いくらでもきりがなく出てくると思われる．あるいはハイパーパラメータと言ってもよいかもしれない．
例
ベイズ線形回帰
TODO
RVM
TODO
EMアルゴリズムが対数尤度を極大化する理由
一般的な議論を行う．以下の尤度関数
$$\begin{align*} p(\mathbf{X} \mid \boldsymbol{\theta})=\sum_{\mathbf{Z}} p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}) \end{align*}$$
を $\boldsymbol{\theta}$ について最大化したい．
$$\begin{align*} \ln p(\mathbf{X} \mid \boldsymbol{\theta})=\mathcal{L}(q, \boldsymbol{\theta})+\mathrm{KL}(q | p) \end{align*}$$
ここで $\mathrm{KL}(q | p)$ はKLダイバージェンスと呼ばれるが常に0以上の値をとる．以下の図は対数尤度を分解した図である．
Eステップでは $q(\mathbf{Z})$ について $\mathcal{L}(q, \boldsymbol{\theta})$ を最大化する．すると下の図のように $\mathrm{KL}(q | p)$ が0になる．
次のMステップでは $\mathcal{L}(q, \boldsymbol{\theta})$ を $\boldsymbol{\theta}$ について最大化する．このとき $\mathrm{KL}(q | p)$ は必ず0以上であるから，図のように $\ln ( \mathbf{X} \mid \boldsymbol{\theta})$ も点線から増加する．
