はじめに
『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでの実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。
この記事は、9.2節の内容です。多次元混合ガウス分布(多変量混合正規分布)の定義の確認と多次元混合ガウス分布に対するEMアルゴリズムによる最尤推定を導出します。
【実装編】
www.anarchive-beta.com
www.anarchive-beta.com
【他の節一覧】
www.anarchive-beta.com
【この節の内容】
9.2.0 混合ガウス分布
$K$個の$D$次元ガウス分布からなる混合ガウス分布を確認する。
$$
p(\mathbf{x})
= \sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\tag{9.7}
$$
ここで、$K$個の平均パラメータを$\boldsymbol{\mu} \equiv \{\boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \cdots, \boldsymbol{\mu}_K\}$、共分散行列パラメータを$\boldsymbol{\Sigma} \equiv \{\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2, \cdots, \boldsymbol{\Sigma}_K\}$とする。$K$個の分布を混合要素と呼ぶ。$\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$はクラスタ$k$の$D$次元ガウス分布であり、平均$\boldsymbol{\mu}_k = (\mu_{k1}, \mu_{k2}, \cdots, \mu_{kD})$、共分散$\boldsymbol{\Sigma}_k = (\sigma_{k11}, \cdots, \sigma_{kDD})$である。
また、$\boldsymbol{\pi} = (\pi_1, \pi_2, \cdots, \pi_K)$を混合係数と呼び、次の条件を満たす。
$$
0 \leq \pi_k \leq 1,\
\sum_{k=1}^K \pi_k = 1
\tag{9.8-9}
$$
つまり、$\boldsymbol{\pi}$は確率の条件を満たしている。
$K$次元の2値確率変数$\mathbf{z} = (z_1, z_2, \cdots, z_K)$を導入して、データ$\mathbf{x}$のクラスタを表す。$\mathbf{z}$の$K$個の要素の内、1つの要素が1で他の要素は0をとる。つまり、$z_k \in \{0, 1\}$であり、$\sum_{k=1}^K z_k = 1$を満たす。このような表現方法を1-of-K符号化法やone-hot表現と呼ぶ。
$z_k = 1$のときクラスタ$k$であることを表す。つまり、クラスタ$k$である確率は
$$
p(z_k = 1)
= \pi_k
$$
である。
($z_k = 1$以外の$z_j = 0\ (j = 1, \cdots, k - 1, k + 1, \cdots, K)$も含めた)$\mathbf{z}$の分布は、パラメータ$\boldsymbol{\pi}$を持つカテゴリ分布で表現できる。
$$
p(\mathbf{z})
= \prod_{k=1}^K
\pi_k^{z_k}
= \mathrm{Cat}(\mathbf{z} | \boldsymbol{\pi})
\tag{9.10}
$$
0乗すると1になるので、クラスタ$k$のとき($z_k = 1$のとき)この式は
$$
\begin{aligned}
p(z_k = 1, z_j = 0)
&= \prod_{k=1}^K
\pi_k^{z_k}
\\
&= \pi_1^{z_1} * \cdots * \pi_k^{z_k} * \cdots * \pi_K^{z_K}
\\
&= \pi_1^0 * \cdots * \pi_k^1 * \cdots * \pi_K^0
\\
&= 1 * \cdots * \pi_k * \cdots * 1
\\
&= \pi_k
\end{aligned}
$$
クラスタ$k$の混合比率$\pi_k$となる($z_k = 1$のとき$z_j = 0$なのは自明なので本来は書かなくてよい)。つまり、$\mathbf{z}$によって対応する確率を取り出せる。
同様に、潜在変数を導入すると、クラスタ$k$のガウス分布とは$z_k =1$を条件とする($z_k = 1$が与えられた下での)ガウス分布なので
$$
p(\mathbf{x} | z_k = 1)
= \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
$$
と書ける。($z_k = 1$が分かっているということは他の項が0なのも分かるので、$\mathbf{z}$が与えられていると言える。)
($z_k = 1$以外の$z_j = 0\ (j \neq k)$も含めた)$\mathbf{z}$が与えられた下でのガウス分布は、次の式で表現できる。
$$
p(\mathbf{x} | \mathbf{z})
= \prod_{k=1}^K
\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)^{z_k}
\tag{9.11}
$$
こちらも$\mathbf{z}$によって対応する分布を取り出せる。
$\mathbf{x}$の分布$p(\mathbf{x})$は、$\mathbf{x}$と$\mathbf{z}$の同時分布
$$
p(\mathbf{x}, \mathbf{z})
= p(\mathbf{x} | \mathbf{z})
p(\mathbf{z})
$$
を$\mathbf{z}$に関して周辺化した分布と言える。
$\mathbf{z}$の全ての組み合わせ($(z_1, z_2, \cdots, z_K)$として$(1, 0, \cdots, 0)$から$(0, 0, \cdots, 1)$まで)の和を$\sum_{\mathbf{z}}$で表現すると、$\mathbf{x}$の周辺分布は次のようになる。
$$
\begin{aligned}
p(\mathbf{x})
&= \sum_{\mathbf{z}}
p(\mathbf{x}, \mathbf{z})
\\
&= \sum_{\mathbf{z}}
p(\mathbf{x} | \mathbf{z})
p(\mathbf{z})
\\
&= \sum_{z}
\prod_{k=1}^K
\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)^{z_k}
\prod_{k=1}^K
\pi_k^{z_k}
\\
&= \sum_{z}
\prod_{k=1}^K \Bigl(
\pi_k
\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\Bigr)^{z_k}
\end{aligned}
$$
$\sum_{z}$と$\prod_{k=1}^K$を展開する。
$$
\begin{aligned}
p(\mathbf{x})
&= \sum_{z} \Bigl\{
\Bigl(
\pi_1 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)
\Bigr)^{z_1}
\Bigl(
\pi_2 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2)
\Bigr)^{z_2}
\cdots
\Bigl(
\pi_K \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_K)
\Bigr)^{z_K}
\Bigr\}
\\
&= \Bigl(
\pi_1 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)
\Bigr)^1
\Bigl(
\pi_2 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2)
\Bigr)^0
\cdots
\Bigl(
\pi_K \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_K)
\Bigr)^0 \\
&\qquad
+ \Bigl(
\pi_1 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)
\Bigr)^0
\Bigl(
\pi_2 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2)
\Bigr)^1
\cdots
\Bigl(
\pi_K \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_K)
\Bigr)^0 \\
&\qquad
\vdots \\
&\qquad
+ \Bigl(
\pi_1 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)
\Bigr)^0
\Bigl(
\pi_2 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2)
\Bigr)^0
\cdots
\Bigl(
\pi_K \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_K)
\Bigr)^1
\\
&= + \pi_1 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)
+ \pi_2 \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2)
+ \cdots
+ \pi_K \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_K, \boldsymbol{\Sigma}_K)
\end{aligned}
$$
1から$K$までの和をまとめる。
$$
p(\mathbf{x})
= \sum_{k=1}^K
\pi_k \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\tag{9.12}
$$
混合ガウス分布の定義式(9.7)となるのが分かる。
最後に、$\mathbf{x}$が与えられた(観測された)下での$\mathbf{z}$の条件付き確率(事後分布)を求める。クラスタ$k$について、尤度関数を$p(\mathbf{x} | z_k = 1)$、事前分布を$p(z_k = 1)$、周辺分布を$p(\mathbf{x})$として、ベイズの定理を用いる。
$$
\begin{align}
p(z_k = 1 | \mathbf{x})
&= \frac{
p(\mathbf{x} | z_k = 1)
p(z_k = 1)
}{
p(\mathbf{x})
}
\\
&= \frac{
\pi_k
\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}{
\sum_{k'=1}^K
\pi_{k'}
\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'})
}
\equiv
\gamma(z_k)
\tag{9.14}
\end{align}
$$
これまでに確認した式をそれぞれ代入すると、混合分布からクラスタ$k$に関して取り出したものを正規化する($k$に関する項を$k = 1, 2, \cdots, K$の総和で割る)形になっている。これを負担率と呼び、$\gamma(z_k)$とおく。この式はEMアルゴリズムにおいて度々登場する。
他のクラスタ$k = 1, \cdots, K$についても同様に求められる。
$\mathbf{z}$の事後分布は、事前分布と同様にカテゴリ分布となる。
$$
p(\mathbf{z} | \mathbf{x})
= \prod_{k=1}^K
\gamma(z_k)^{z_k}
= \mathrm{Cat}(\mathbf{z} | \gamma(\mathbf{z}))
$$
ここで、パラメータを$\gamma(\mathbf{z}) = (\gamma(z_{1}), \gamma(z_1), \cdots, \gamma(z_K))$とおき、次の条件を満たす。
$$
0 \leq \gamma(z_k) \leq 1,\
\sum_{k=1}^K \gamma(z_k) = 1
$$
9.2.1 最尤推定
$N$個の観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$がそれぞれ独立に生成されると仮定すると、尤度関数は次のように分解できる。
$$
\begin{aligned}
p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
&= \prod_{n=1}^N
p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
\\
&= \prod_{n=1}^N
\sum_{k=1}^K
\pi_k
p(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\end{aligned}
$$
対数をとった尤度を対数尤度と呼ぶ。
$$
\begin{align}
\ln p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
&= \sum_{n=1}^N
\ln p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
\\
&= \sum_{n=1}^N
\ln \left\{
\sum_{k=1}^K
\pi_k
p(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\right\}
\equiv
L
\tag{9.14}
\end{align}
$$
次項では、対数尤度$L$を最大化する各パラメータを求めていく。
9.2.2 混合ガウス分布のEMアルゴリズム
EMアルゴリズムを用いて、尤度関数が最大となるパラメータ$\boldsymbol{\mu},\ \boldsymbol{\Sigma},\ \boldsymbol{\pi}$を求める。尤度を各パラメータで偏微分して0となる(極大値となる)パラメータを最尤解(最尤推定量)と呼ぶ。ただし、計算を簡単にするため対数をとった尤度関数を用いる。対数関数は単調増加するため結果は一致する。
・平均の最尤解の導出
まずは、対数尤度を最大化する平均$\boldsymbol{\mu}$を求めていく。
対数尤度$L$を$\boldsymbol{\mu}_k$に関して微分する。ただし、式が複雑なため1つずつ確認していく。
$$
\begin{align}
\frac{\partial L}{\partial \boldsymbol{\mu}_k}
&= \frac{\partial}{\partial \boldsymbol{\mu}_k}
\ln p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
\\
&= \sum_{n=1}^N
\frac{\partial}{\partial \boldsymbol{\mu}_k}
\ln p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
\tag{1}
\end{align}
$$
$N$個のデータは独立しているため、対数尤度の微分は、$N$個の対数をとった混合ガウス分布の微分の和になる。
$n$番目のデータに関する対数をとった混合ガウス分布を微分する。
$$
\frac{\partial}{\partial \boldsymbol{\mu}_k}
\ln p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
= \frac{1}{p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})}
\frac{\partial}{\partial \boldsymbol{\mu}_k}
p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
$$
---------- ここから補足 ----------
対数をとった関数$\ln f(x)$の微分を考える。対数関数$\ln$を$g(\cdot)$で表すと、$\ln f(x)$は合成関数$g(f(x))$で表現できる。よって$\ln f(x)$の微分は、合成関数の微分$\{g(f(x))\}' = g'(f(x)) f'(x)$より、「$f(x)$による$g(f(x))$の微分$g'(f(x))$」と「$x$による$f(x)$の微分$f'(x)$」の積になる。さらに$g'(f(x))$は、自然対数の微分$(\ln x)' = \frac{1}{x}$より、$\frac{1}{f(x)}$になる。したがって$\ln f(x)$の微分は、$\frac{f'(x)}{f(x)}$である。つまり、次のように書ける。
$$
\frac{\partial \ln f(x)}{\partial x}
= \frac{1}{f(x)}
\frac{\partial f(x)}{\partial x}
$$
確率密度関数$p(\cdot)$を$f(\mu_k)$とし、対数関数(式全体)$\ln p(\cdot)$を$g(f(\mu_k))$として式を変形する。
---------- ここまで ----------
(対数をとっていない)混合ガウス分布(9.7)を微分する。
$$
\begin{aligned}
\frac{\partial}{\partial \boldsymbol{\mu}_k}
p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
&= \frac{\partial}{\partial \boldsymbol{\mu}_k}
\sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\\
&= \pi_k
\frac{\partial}{\partial \boldsymbol{\mu}_k}
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\end{aligned}
$$
---------- ここから補足 ----------
$\boldsymbol{\mu}_k$に関して微分するので、$\boldsymbol{\mu}_k$を含まない$1, \cdots, k - 1, k + 1, \cdots, K$の項は次のように0になり消える。
$$
\begin{aligned}
\frac{\partial}{\partial x_k}
\sum_{k=1}^K f(x_k)
&= \frac{\partial f(x_1)}{\partial x_k} + \cdots
+ \frac{\partial f(x_k)}{\partial x_k}
+ \cdots + \frac{\partial f(x_K)}{\partial x_k}
\\
&= 0 + \cdots
+ \frac{\partial f(x_k)}{\partial x_k}
+ \cdots + 0
\end{aligned}
$$
---------- ここまで ----------
クラスタ$k$の多次元ガウス分布の微分を対数微分法を用いて求める。
$$
\frac{\partial}{\partial \boldsymbol{\mu}_k}
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
= \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\frac{\partial}{\partial \boldsymbol{\mu}_k}
\ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
$$
---------- ここから補足 ----------
先ほど確認した対数をとった関数の微分$\{\ln f(x)\}' = \frac{f'(x)}{f(x)}$について、両辺に$f(x)$を掛けると、元の関数の微分$f'(x)$は「元の関数$f(x)$」と「対数をとった関数の微分$\{\ln f(x)\}'$」の積によって求められることが分かる。
$$
\frac{\partial f(x)}{\partial x}
= f(x)
\frac{\partial \ln f(x)}{\partial x}
$$
これを対数微分法と呼び、$f(x)$が複雑な関数で$\ln f(x)$の微分の方が簡単にできるとき用いられる。
---------- ここまで ----------
多次元ガウス分布の微分の具体的な式を求める。クラスタ$k$の多次元ガウス分布は次の式である。
$$
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
= \frac{1}{(2 \pi)^{\frac{D}{2}}}
\frac{1}{|\boldsymbol{\Sigma}_k|^{\frac{1}{2}}}
\exp \left\{
- \frac{1}{2}
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
\right\}
$$
この式の対数をとり$\boldsymbol{\mu}_k$に関して整理する。
$$
\begin{aligned}
\ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
&= - \frac{1}{2}
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
+ \mathrm{const.}
\\
&= - \frac{1}{2} \Bigl(
\mathbf{x}_n^{\top} \boldsymbol{\Sigma}_k^{-1} \mathbf{x}_n
- 2 \boldsymbol{\mu}_k^{\top} \boldsymbol{\Sigma}_k^{-1} \mathbf{x}_n
+ \boldsymbol{\mu}_k^{\top} \boldsymbol{\Sigma}_k^{-1} \boldsymbol{\mu}_k
\Bigr)
+ \mathrm{const.}
\\
&= \boldsymbol{\mu}_k^{\top} \boldsymbol{\Sigma}_k^{-1} \mathbf{x}_n
- \frac{1}{2}
\boldsymbol{\mu}_k^{\top} \boldsymbol{\Sigma}_k^{-1} \boldsymbol{\mu}_k
+ \mathrm{const.}
\end{aligned}
$$
$\boldsymbol{\mu}_k$に影響しない項を$\mathrm{const.}$にまとめた。$\mathrm{const.}$は$\boldsymbol{\mu}_k$と無関係なので微分時に0になる。
この式を$\boldsymbol{\mu}_k$に関して微分する。
$$
\begin{aligned}
\frac{\partial}{\partial \boldsymbol{\mu}_k}
\ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
&= \boldsymbol{\Sigma}_k^{-1} \mathbf{x}_n
- \frac{1}{2} \Bigl\{
\boldsymbol{\Sigma}_k^{-1} + (\boldsymbol{\Sigma}_k^{-1})^{\top}
\Bigr\}
\boldsymbol{\mu}_k
\\
&= \boldsymbol{\Sigma}_k^{-1}
\mathbf{x}_n
- \boldsymbol{\Sigma}_k^{-1}
\boldsymbol{\mu}_k
\\
&= \boldsymbol{\Sigma}_k^{-1}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
\end{aligned}
$$
---------- ここから補足 ----------
前の項は、$\boldsymbol{\Sigma}_k^{-1} \mathbf{x}_n$を1つのベクトル$\mathbf{a}$として付録Cの式(C.19)、$\frac{\partial \mathbf{x}^{\top} \mathbf{a}}{\partial \mathbf{x}} = \mathbf{a}$の変形を行う。
後の項は、「The Matrix Cookbook」の式(81)、$\frac{\partial \mathbf{x} \mathbf{B} \mathbf{x}}{\partial \mathbf{x}} = (\mathbf{B} + \mathbf{B}^{\top}) \mathbf{x}$の変形を行う。また、$\boldsymbol{\Sigma}_k^{-1}$は対称行列なので$\boldsymbol{\Sigma}_k^{-1} = (\boldsymbol{\Sigma}_k^{-1})^{\top}$である。
ちなみに、この式全体の微分も式(85)として載っている。
---------- ここまで ----------
これまで求めた各項の微分を式(1)に代入すると、対数尤度$L$の$\boldsymbol{\mu}_k$に関する微分が得られる。
$$
\begin{aligned}
\frac{\partial L}{\partial \boldsymbol{\mu}_k}
&= \sum_{n=1}^N
\frac{1}{p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})}
\pi_k
p(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\boldsymbol{\Sigma}_k^{-1}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
\\
&= \sum_{n=1}^N
\frac{
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}{
\sum_{k'=1}^K
\pi_{k'}
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'})
}
\boldsymbol{\Sigma}_k^{-1}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
\\
&= \boldsymbol{\Sigma}_k^{-1}
\sum_{n=1}^N
\gamma(z_{nk})
(\mathbf{x}_n - \boldsymbol{\mu}_k)
\end{aligned}
$$
分数の部分が負担率の定義式(9.13)なので$\gamma(z_{nk})$に置き換えた。
対数尤度の微分$\frac{\partial L}{\partial \boldsymbol{\mu}_k}$を0とおき、$\boldsymbol{\mu}_k$に関して解く。
$$
\begin{align}
&\frac{\partial L}{\partial \boldsymbol{\mu}_k}
= \boldsymbol{\Sigma}_k^{-1}
\sum_{n=1}^N
\gamma(z_{nk})
(\mathbf{x}_n - \boldsymbol{\mu}_k)
= 0
\tag{9.16}\\
&\Leftrightarrow
\sum_{n=1}^N
\gamma(z_{nk})
\mathbf{x}_n
- \boldsymbol{\mu}_k
\sum_{n=1}^N
\gamma(z_{nk})
= 0
\\
&\Leftrightarrow
\boldsymbol{\mu}_k
\sum_{n=1}^N
\gamma(z_{nk})
= \sum_{n=1}^N
\gamma(z_{nk}) \mathbf{x}_n
\end{align}
$$
1行目では、両辺に左から$\boldsymbol{\Sigma}_k$を掛けて$\boldsymbol{\Sigma}_k^{-1}$を打ち消した。ここで
$$
N_k = \sum_{n=1}^N
\gamma(z_{nk})
\tag{9.18}
$$
とおく。$\gamma(z_{nk})$は$n$番目のデータがクラスタ$k$となる確率を表す。よって、$N_k$はクラスタ$k$となるデータ数の期待値と言える。
したがって、$\boldsymbol{\mu}_k$の最尤解は次の式となる。
$$
\boldsymbol{\mu}_k
= \frac{1}{N_k}
\sum_{n=1}^N
\gamma(z_{nk}) \mathbf{x}_n
\tag{9.17}
$$
各データの負担率による期待値と言える($N$個分のデータの和を求めているが、実質的には$N_k$個分の和である)。
他のクラスタ$k = 1, \cdots, K$についても同様に求められる。
・共分散行列の最尤解の導出
同様に、対数尤度を最大化する共分散行列$\boldsymbol{\Sigma}$を求めていく。
対数尤度$L$を$\boldsymbol{\Sigma}_k$に関して微分する。
$$
\begin{align}
\frac{\partial L}{\partial \boldsymbol{\Sigma}_k}
&= \frac{\partial}{\partial \boldsymbol{\Sigma}_k}
\ln p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
\\
&= \sum_{n=1}^N
\frac{\partial}{\partial \boldsymbol{\Sigma}_k}
\ln p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
\\
&= \sum_{n=1}^N
\frac{1}{p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})}
\frac{\partial}{\partial \boldsymbol{\Sigma}_k}
p(\mathbf{x}_n | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
\\
&= \sum_{n=1}^N
\frac{
\pi_k
}{
\sum_{k'=1}^K
\pi_{k'}
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'})
}
\frac{\partial}{\partial \boldsymbol{\Sigma}_k}
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\\
&= \sum_{n=1}^N
\frac{
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}{
\sum_{k'=1}^K
\pi_{k'}
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_{k'}, \boldsymbol{\Sigma}_{k'})
}
\frac{\partial}{\partial \boldsymbol{\Sigma}_k}
\ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\tag{2}
\end{align}
$$
ここまでは平均パラメータのときと同じ手順である。
対数をとったクラスタ$k$の多次元ガウス分布の微分を考える。今度は$\boldsymbol{\Sigma}_k$に関して式を整理する。
$$
\ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
= - \frac{1}{2}
\ln |\boldsymbol{\Sigma}_k|
- \frac{1}{2}
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
+ \mathrm{const.}
$$
$\boldsymbol{\Sigma}_k$と無関係な項を$\mathrm{const.}$にまとめる。
この式を$\boldsymbol{\Sigma}_k$に関して微分する。
$$
\begin{aligned}
\frac{\partial}{\partial \boldsymbol{\Sigma}_k}
\ln \mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
&= - \frac{1}{2}
(\boldsymbol{\Sigma}_k^{-1})^{\top}
+ \frac{1}{2}
(\boldsymbol{\Sigma}_k^{-1})^{\top}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
(\boldsymbol{\Sigma}_k^{-1})^{\top}
\\
&= - \frac{1}{2}
\boldsymbol{\Sigma}_k^{-1}
+ \frac{1}{2}
\boldsymbol{\Sigma}_k^{-1}
(\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
\\
&= - \frac{1}{2}
\boldsymbol{\Sigma}_k^{-1} \Bigl\{
\mathbf{I}
- (\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
\Bigr\}
\end{aligned}
$$
---------- ここから補足 ----------
前の項は、付録Cの式(C.28)、$\frac{\partial}{\partial \mathbf{A}} |\mathbf{A}| = (\mathbf{A}^{-1})^{\top}$の変形を行う。また、$\boldsymbol{\Sigma}_k^{-1}$は対称行列なので、$\boldsymbol{\Sigma}_k^{-1} = (\boldsymbol{\Sigma}_k^{-1})^{\top}$である。
後の項は、「The Matrix Cookbook」の式(61)、$\frac{\partial \mathbf{a}^{\top} \mathbf{X}^{-1} \mathbf{b}}{\partial \mathbf{X}} = - (\mathbf{X}^{-1})^{\top} \mathbf{a} \mathbf{b}^{\top} (\mathbf{X}^{-1})^{\top}$の変形を行う。
---------- ここまで ----------
この式を対数尤度の微分(2)に代入すると、対数尤度$L$の$\boldsymbol{\Sigma}_k$に関する微分が得られる。
$$
\begin{aligned}
\frac{\partial L}{\partial \boldsymbol{\Sigma}_k}
&= \sum_{n=1}^N
\frac{
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}{
\sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}
\left(
- \frac{1}{2}
\boldsymbol{\Sigma}_k^{-1}
\right) \Bigl\{
\mathbf{I}
- (\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
\Bigr\}
\\
&= - \frac{1}{2}
\boldsymbol{\Sigma}_k^{-1}
\sum_{n=1}^N
\gamma(z_{nk}) \Bigl\{
\mathbf{I}
- (\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
\Bigr\}
\end{aligned}
$$
負担率の定義式(9.13)より、$\gamma(z_{nk})$に置き換えた。
対数尤度の微分$\frac{\partial L}{\partial \boldsymbol{\Sigma}_k}$を0とおき、$\boldsymbol{\Sigma}_k$に関して解く。
$$
\begin{align}
&\frac{\partial L}{\partial \boldsymbol{\Sigma}_k}
= - \frac{1}{2}
\boldsymbol{\Sigma}_k^{-1}
\sum_{n=1}^N
\gamma(z_{nk}) \Bigl\{
\mathbf{I}
- (\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
\Bigr\}
= 0
\\
&\Leftrightarrow
\sum_{n=1}^N
\gamma(z_{nk})
\mathbf{I}
- \sum_{n=1}^N
\gamma(z_{nk})
(\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
= 0
\\
&\Leftrightarrow
N_k
\mathbf{I}
= \sum_{n=1}^N
\gamma(z_{nk})
(\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\boldsymbol{\Sigma}_k^{-1}
\\
&\Leftrightarrow
\boldsymbol{\Sigma}_k
= \frac{1}{N_k}
\sum_{n=1}^N
\gamma(z_{nk})
(\mathbf{x}_n - \boldsymbol{\mu}_k)
(\mathbf{x}_n - \boldsymbol{\mu}_k)^{\top}
\tag{9.19}
\end{align}
$$
1行目では、両辺に左から$- 2 \boldsymbol{\Sigma}_k$を掛けて$\boldsymbol{\Sigma}_k^{-1}$を打ち消し、また波括弧を展開して$\boldsymbol{\Sigma}_k$に関して式を整理する。最後の式変形では、両辺に右から$\boldsymbol{\Sigma}_k$を掛けて最尤解が得られる。
偏差の2乗の負担率による期待値と言える。
他のクラスタ$k = 1, \cdots, K$についても同様に求められる。
・混合係数の最尤解の導出
最後に、対数尤度を最大化する混合係数$\boldsymbol{\pi}$を求めていく。
$\boldsymbol{\pi}$には総和が1であるという制約条件(9.9)がある。そこで、ラグランジュ未定係数法を用いて、対数尤度に制約条件を含めた式を立てる。
$$
\begin{align}
F(\boldsymbol{\pi})
&= \ln p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})
+ \lambda \left(
\sum_{k=1}^K \pi_k - 1
\right)
\tag{9.20}\\
&= \sum_{n=1}^N
\ln \left(
\sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\right)
+ \lambda \left(
\sum_{k=1}^K \pi_k - 1
\right)
\end{align}
$$
立てた式を$F(\boldsymbol{\pi})$とする。
これまでのように、$F(\boldsymbol{\pi})$を$\pi_k$に関して微分する。
$$
\begin{aligned}
\frac{\partial F(\boldsymbol{\pi})}{\partial \pi_k}
&= \sum_{n=1}^N
\frac{\partial}{\partial \pi_k}
\ln \left(
\sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
\right)
+ \frac{\partial}{\partial \pi_k}
\lambda \left(
\sum_{k=1}^K \pi_k - 1
\right)
\\
&= \sum_{n=1}^N
\frac{
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}{
\sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}
+ \lambda
\end{aligned}
$$
---------- ここから補足 ----------
前の因子の微分は、これまでのように括弧全体を$f(\pi_k)$として対数をとった関数の微分を行う。$f(\pi_k)$の微分は、$\pi_k$を含まないガウス分布$\mathcal{N}(\cdot)$を定数として扱う。
後の因子の微分は、次のように$\pi_k$に影響する項以外は消える。
$$
\begin{aligned}
\frac{\partial}{\partial \pi_k}
\lambda \left(
\sum_{k=1}^K \pi_k - 1
\right)
&= \frac{\partial}{\partial \pi_k} \Bigl(
\lambda \pi_1 + \cdots + \lambda \pi_k + \cdots + \lambda \pi_K - 1
\Bigr)
\\
&= 0 + \cdots + \lambda + 0 + 0
\end{aligned}
$$
---------- ここまで ----------
$F(\pi_k)$の微分を0とおき、$\pi_k$に関して解く。
$$
\begin{align}
&\frac{\partial F(\boldsymbol{\pi})}{\partial \pi_k}
= \sum_{n=1}^N
\frac{
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}{
\sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}
+ \lambda
= 0
\tag{9.21}\\
&\Leftrightarrow
\sum_{n=1}^N
\frac{
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}{
\sum_{k=1}^K
\pi_k
\mathcal{N}(\mathbf{x}_n | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)
}
+ \lambda \pi_k
= 0
\end{align}
$$
式(9.21)の両辺に$\pi_k$を掛けると、前の因子が負担率の定義式(9.13)になる。$\gamma(z_{nk})$に置き換える。
$$
\begin{aligned}
\sum_{n=1}^N
\gamma(z_{nk})
+ \lambda \pi_k
&= 0
\\
\Leftrightarrow
\lambda \pi_k
&= - \sum_{n=1}^N
\gamma(z_{nk})
\\
&= - N_k
\end{aligned}
$$
両辺で$k$について和をとる。
$$
\begin{aligned}
\lambda
\sum_{k=1}^K \pi_k
&= - \sum_{k=1}^K
N_k
\\
\Leftrightarrow
\lambda
&= - N
\end{aligned}
$$
混合比率の制約条件より$\sum_{k=1}^K \pi_k = 1$、また負担率の定義より$\sum_{n=1}^N \sum_{k=1}^K \gamma(z_{nk}) = N$である。
この式を1つ上の式に代入すると、$\pi_k$の最尤解が得られる。
$$
\begin{align}
- N \pi_k
&= - N_k
\\
\Leftrightarrow
\pi_k
&= \frac{N_k}{N}
\tag{9.22}
\end{align}
$$
クラスタ$k$となる確率$\pi_k$は、各データに割り当てられたクラスタ$k$となる確率の総和$N_k$が全データに占める割合によって決まることが分かる。
他のクラスタ$k = 1, \cdots, K$についても同様に求められる。
以上で、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Sigma},\ \boldsymbol{\pi}$の最尤解が得られた。Eステップでは、前ステップ(あるいは初期値)の$\boldsymbol{\mu},\ \boldsymbol{\Sigma},\ \boldsymbol{\pi}$を用いて、負担率$\gamma(z_{nk})$を更新する。ちなみに、$\mathbf{z}$の事後分布はカテゴリ分布なので、$\mathbf{z}$の期待値は負担率$\mathbb{E}_{p(\mathbf{z} | \mathbf{x})}[\mathbf{z}] = \gamma(\mathbf{z})$である。Mステップでは、更新した$\gamma(z_{nk})$を用いて、パラメータ$\boldsymbol{\mu},\ \boldsymbol{\Sigma},\ \boldsymbol{\pi}$を更新する。更新したパラメータは尤度関数を増加させる。これを繰り返すことで尤度関数を最大値に近付けていく。
参考文献
- C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.
- K.B.Petersen, M.S.Pedersen. The Matrix Cookbook. Technical University of Denmark, 2012.
おわりに
ベイズ推論による機械学習入門を読み終わったのでこの本も少しずつ進めます。
【関連する内容】
変分推論はこちらでやりました。
www.anarchive-beta.com
