混合正規分布(Mixtures of Gaussians)の尤度関数とEMアルゴリズムの導出

混合正規分布(Mixtures of Gaussians)は多峰の確率分布の表現にあたって、複数の正規分布を確率的に混合して表す分布です。当記事では混合正規分布の尤度関数を確認し、尤度最大化にあたって用いるEMアルゴリズムの導出について確認を行いました。

「パターン認識と機械学習」の$2.3.9$節と$9.2$節の「混合正規分布(Mixtures of Gaussians)」を参考に作成を行いました。

混合正規分布

混合正規分布の式定義

$$
\large
\begin{align}
p(x) = \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x|\mu_k,\Sigma_k) \quad (1)
\end{align}
$$

ここで$\pi_{k}$は混合係数(mixing coefficients)であり、分布を混ぜ合わせる割合であることから$\pi_{k}$は下記の制約が成立する必要がある。
$$
\large
\begin{align}
\sum_{k=1}^{K} \pi_{k} &= 1 \\
0 \leq \pi_{k} \leq 1
\end{align}
$$

このとき、$(1)$式は下記のように表せることも抑えておくとよい。
$$
\large
\begin{align}
p(x) &= \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x|\mu_k,\Sigma_k) \\
&= \sum_{k=1}^{K} p(k)p(x|k) \quad (2)
\end{align}
$$

負担率の定義

$(2)$式の表記では$\displaystyle p(k)=\pi_{k}, p(x|k)=\mathcal{N}(x|\mu,\Sigma)$のように対応させているが、ここで$k$の事後分布$p(k|x)$を負担率$\gamma_{k} \equiv p(k|x)$のように表すと、パラメータ推定の表記を簡易化できる。

負担率$\gamma_{k}(x) \equiv p(k|x)$の式は下記のように考えることができる。
$$
\large
\begin{align}
\gamma_{k}(x) & \equiv p(k|x) \\
&= \frac{p(k)p(x|k)}{\displaystyle \sum_{l=1}^{K} p(l)p(x|l)} \\
&= \frac{\pi_{k} \mathcal{N}(x|\mu,\Sigma)}{\displaystyle \sum_{l=1}^{K} \pi_{l} \mathcal{N}(x|\mu_l,\Sigma_l)} \quad (3)
\end{align}
$$

混合正規分布の尤度関数

混合正規分布の尤度関数を$\mathcal{L}(\pi,\mu,\Sigma)$とおくと、$\mathcal{L}(\pi,\mu,\Sigma)$は観測値$x_1,…,x_N$に関する同時確率密度関数に一致するので下記のように考えることができる。
$$
\large
\begin{align}
\mathcal{L}(\pi,\mu,\Sigma) &= \prod_{n=1}^{N} p(x_n) \\
&= \prod_{n=1}^{N} \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k)
\end{align}
$$

上記より、対数尤度$\log{\mathcal{L}(\pi,\mu,\Sigma)}$を下記のように考えられる。
$$
\large
\begin{align}
\log{\mathcal{L}(\pi,\mu,\Sigma)} &= \log{ \prod_{n=1}^{N} \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k) } \\
&= \log{ \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_1|\mu_k,\Sigma_k) \times … \times \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k) } \\
&= \sum_{n=1}^{N} \log{ \left[ \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k) \right] } \quad (4)
\end{align}
$$

EMアルゴリズム

潜在変数$z$の導入

$(2)$式を下記のように定義する$1$-of-$K$表現の潜在変数$z$を用いて書き換えることを考える。
$$
\large
\begin{align}
z &= \left(\begin{array}{c} z_{1} \\ \vdots \\ z_{K} \end{array} \right) \\
z_k & \in \{ 0,1 \} \\
\sum_{k=1}^{K} z_k &= 1
\end{align}
$$

数式だとやや難しく見えるが、上記は「$z$は$z_1$から$z_K$のどれか$1$つが$1$で他が全て$0$の$K$次元ベクトルである」と解釈できるので解釈自体はそれほど難しくはない。ここで$p(z_k=1)=\pi_k$であると考えると、$p(z)$は下記のように考えることができる。
$$
\large
\begin{align}
p(z) &= p(z_1,z_2,…,z_K) \\
&= \prod_{k=1}^{K} p(z_k) \\
&= \prod_{k=1}^{K} \pi_k^{z_k} \quad (5)
\end{align}
$$

同時確率密度関数$p(x,z)$を考えるとき、$(5)$式は$z$に関する周辺分布の確率密度関数であると考えることができる。また、同様に条件付き確率$p(x|z)$も下記のように考えられる。
$$
\large
\begin{align}
p(x|z) = \prod_{k=1}^{K} \mathcal{N}(x|\mu_k,\Sigma_k)^{z_k}
\end{align}
$$

よって$p(x,z)$に対する$x$の周辺分布$p(x)$は下記のように表すことができる。
$$
\large
\begin{align}
p(x) &= \sum_{z} p(z)p(x|z) \\
&= \sum_{k=1}^{K} \pi_k \mathcal{N}(x|\mu_k,\Sigma_k) \quad (1)’
\end{align}
$$

上記の詳しい計算は章末演習$9$.$3$のように考えると良い。$(1)’$式より対数尤度は下記のように$(4)$式と同様に考えられる。
$$
\large
\begin{align}
\log{\mathcal{L}(\pi,\mu,\Sigma)} = \sum_{n=1}^{N} \log{ \left[ \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k) \right] } \quad (6)
\end{align}
$$

また、負担率$\gamma(z_k)$を下記のように定義する。
$$
\large
\begin{align}
\gamma(z_k) & \equiv p(z_k=1|x) \\
&= \frac{\pi_{k} \mathcal{N}(x|\mu,\Sigma)}{\displaystyle \sum_{l=1}^{K} \pi_{l} \mathcal{N}(x|\mu_l,\Sigma_l)} \quad (7)
\end{align}
$$

EMアルゴリズムの導出

$(6)$式を$\mu_k$に関して偏微分を行うと下記が得られる。
$$
\large
\begin{align}
\frac{\partial}{\partial \mu_k} & \log{\mathcal{L}(\pi,\mu,\Sigma)} = \frac{\partial}{\partial \mu_k} \sum_{n=1}^{N} \log{ \left[ \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k) \right] } \\
&= \sum_{n=1}^{N} \frac{\partial}{\partial \mu_k} \log{ \left[ \sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k) \right] } \\
&= \sum_{n=1}^{N} \frac{\displaystyle \frac{\partial}{\partial \mu_k} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k)}{\sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k)} \\
&= \sum_{n=1}^{N} \frac{\mathcal{N}(x_n|\mu_k,\Sigma_k)}{\sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k)} \frac{\partial}{\partial \mu_k} \left( -\frac{1}{2}(x_n-\mu_k)^{\mathrm{T}} \Sigma_k^{-1} (x_n-\mu_k) \right) \\
&= -\sum_{n=1}^{N} \frac{\mathcal{N}(x_n|\mu_k,\Sigma_k)}{\sum_{k=1}^{K} \pi_{k} \mathcal{N}(x_n|\mu_k,\Sigma_k)} \Sigma_k^{-1} (x_n-\mu_k) \\
&= -\sum_{n=1}^{N} \gamma(z_{nk}) \Sigma_k^{-1} (x_n-\mu_k)
\end{align}
$$

ここで「上記$=0$」を$\mu_k$に関して解くと下記が得られる。
$$
\large
\begin{align}
\frac{\partial}{\partial \mu_k} \log{\mathcal{L}(\pi,\mu,\Sigma)} &= 0 \\
\sum_{n=1}^{N} \gamma(z_{nk}) \Sigma_k^{-1} (x_n-\mu_k) &= 0 \\
\sum_{n=1}^{N} \gamma(z_{nk}) (x_n-\mu_k) &= 0 \\
\sum_{n=1}^{N} \gamma(z_{nk}) \mu_k &= \sum_{n=1}^{N} \gamma(z_{nk}) x_n \\
N_k \mu_k &= \sum_{n=1}^{N} \gamma(z_{nk}) x_n \\
\mu_k &= \frac{1}{N_k} \sum_{n=1}^{N} \gamma(z_{nk}) x_n \quad (8)
\end{align}
$$

上記を計算するにあたって、$\displaystyle N_{k} = \sum_{n=1}^{N} \gamma(z_{nk})$のようにおいたが、「$N_k$は負担率$\gamma(z_{nk})$が表す各サンプルの事後確率$p(z_{nk}|x_n)$の$N$個のサンプル分の和」と解釈できる。

$\Sigma_k, \pi_{k}$に関しても同様に考えると下記が導出できる。
$$
\large
\begin{align}
\Sigma_k &= \frac{1}{N_k} \sum_{n=1}^{N} \gamma(z_{nk}) (x_n-\mu_k)(x_n-\mu_k)^{\mathrm{T}} \quad (9) \\
\pi_{k} &= \frac{N_k}{N} \quad (10)
\end{align}
$$
ただし、$\pi_{k}$の計算にあたっては制約条件$\displaystyle \sum_{k=1}^{K} \pi_{k} = 1$を考慮するにあたってラグランジュの未定乗数法を用いる。

$(8), (9), (10)$式には$\gamma(z_{nk})$があることでこの式だけでは最適解を考えることができないが、$(7)$式を元に繰り返しのスキームを用いることで最適解を考えることができる。ここで$(7)$式をEステップ、$(8), (9), (10)$式をMステップと考えるとこのスキームはEMアルゴリズムに合致する。このようにEMアルゴリズムでは二つのステップを交互に繰り返すことで最適解の計算を行う。

EMアルゴリズムの手順まとめ

以下の手順に沿ってEMアルゴリズムを元にパラメータ推定を行うことができる。

$1$. $\pi_k, \mu_k,\Sigma_k$に関してそれぞれ初期値を設定する。

$2$. $\pi_k, \mu_k,\Sigma_k$を固定し、$(7)$式に沿って$\gamma(z_{nk})$の計算を行う。$\leftarrow$Eステップ

$3$. $\gamma(z_{nk})$を固定し、$(8), (9), (10)$式に沿って$\pi_k, \mu_k,\Sigma_k$の計算を行う。$\leftarrow$Mステップ

$4$. $(6)$式に沿って対数尤度を計算し、収束していなければ$2$に戻り処理を繰り返す。

参考

・連続型確率分布
https://www.hello-statisticians.com/explain-terms-cat/probdist2.html

「混合正規分布(Mixtures of Gaussians)の尤度関数とEMアルゴリズムの導出」への1件の返信

コメントは受け付けていません。