からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

9.3.3:混合ベルヌーイ分布のEMアルゴリズム【PRMLのノート】

はじめに

 『パターン認識と機械学習』の独学時のまとめです。一連の記事は「数式の行間埋め」または「R・Pythonでの実装」からアルゴリズムの理解を補助することを目的としています。本とあわせて読んでください。

 この記事は、9.3.3項の内容です。混合ベルヌーイ分布の定義の確認と、混合ベルヌーイ分布に対するEMアルゴリズムによる最尤推定を導出します。

【実装編】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

9.3.3 混合ベルヌーイ分布

 この項では、$N$個の観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$が$K$個のパラメータ$\boldsymbol{\mu} = \{\boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \cdots, \boldsymbol{\mu}_K\}$を持つ混合ベルヌーイ分布に従うとする。各データは$D$次元ベクトル$\mathbf{x}_n = (x_{n1}, x_{n2}, \cdots, x_{nD})^{\top}$であり、1つのクラスタ$k = 1, 2, \cdots, K$が割り当てられているとする。また$\mathbf{x}_n$は、割り当てられたクラスタのパラメータ$\boldsymbol{\mu}_k = (\mu_{k1}, \mu_{k2}, \cdots, \mu_{kD})^{\top}$を持つ$D$個のベルヌーイ分布によって独立に生成されるとする。

・混合ベルヌーイ分布の定義

 まずは、混合ベルヌーイ分布について確認する。

 各データ$\mathbf{x}_n$は、次の混合ベルヌーイ分布によって独立に生成されると仮定する。

$$ \begin{align} p(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) &= \sum_{k=1}^K \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \tag{9.47}\\ &= \sum_{k=1}^K \pi_k \prod_{i=1}^D \mathrm{Bern}(\mathbf{x}_n | \boldsymbol{\mu}_k) \\ &= \sum_{k=1}^K \pi_k \prod_{i=1}^D \mu_{ki}^{x_{ni}} (1 - \mu_{ki})^{1 - x_{ni}} \end{align} $$

 ここで、混合係数を$\boldsymbol{\pi} = (\pi_1, \pi_2, \cdots, \pi_K)^{\top}$、$0 \leq \pi \leq 1,\ \sum_{k=1}^K \pi_k = 1$として、$\mathbf{x}_n$がクラスタ$k$となる確率を$\pi_k$で表す。

 観測データ$\mathbf{X}$の各要素$x_{ni}$は、0か1の値をとる2値変数$x_{ni} \in \{0, 1\}$であり、$\mathbf{x}_n$にクラスタ$k$が割り当てられたとき、パラメータ$\mu_{ki}$を持つベルヌーイ分布に従うとする。

$$ p(x_{ni} | \mu_{ki}) = \mathrm{Bern}(x_{ni} | \mu_{ki}) = \mu_{ki}^{x_{ni}} (1 - \mu_{ki})^{1 - x_{ni}} \tag{2.2} $$

 また、$\mathbf{x}_n$にクラスタ$k$が割り当てられたとき、$\mathbf{x}_n$は$\boldsymbol{\mu}_k$をパラメータとする$D$個のベルヌーイ分布に従うとする。

$$ p(\mathbf{x}_n | \boldsymbol{\mu}_k) = \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) = \prod_{i=1}^D \mu_{ki}^{x_{ni}} (1 - \mu_{ki})^{1 - x_{ni}} \tag{9.44} $$


 混合ガウス分布(9.2節)と同様に、各データのクラスタを示す潜在変数$\mathbf{Z} = \{\mathbf{z}_1, \mathbf{z}_2, \cdots, \mathbf{z}_N\}$、$\mathbf{z}_n = (z_{n1}, z_{n2}, \cdots, z_{nK})^{\top}$を導入する。$z_{nk} \in \{0, 1\},\ \sum_{k=1}^K z_{nk} = 1$であり、$z_{nk} = 1$のとき$\mathbf{x}_n$のクラスタが$k$であることを表す。

 $\mathbf{x}_n$の潜在変数$\mathbf{z}_n$が与えられたとき、$\mathbf{x}_n$の条件付き分布は次のように表せる。

$$ \begin{align} p(\mathbf{x}_n | \mathbf{z}_n, \boldsymbol{\mu}) &= \prod_{k=1}^K p(\mathbf{x}_n | \boldsymbol{\mu}_k)^{z_{nk}} \tag{9.52}\\ &= \prod_{k=1}^K \left\{ \prod_{i=1}^D \mathrm{Bern}(\mathbf{x}_n | \boldsymbol{\mu}_k) \right\}^{z_{nk}} \end{align} $$

 $a^0 = 1$なので、$z_{nk} = 1$のときこの分布は

$$ p(\mathbf{x}_n | z_{nk} = 1, \boldsymbol{\mu}_k) = p(\mathbf{x}_n | \boldsymbol{\mu}_k) \tag{9.44} $$

となる。

 $\mathbf{z}_n$は、混合係数$\boldsymbol{\pi}$をパラメータとするカテゴリ分布に従う。

$$ p(\mathbf{z}_n | \boldsymbol{\pi}) = \mathrm{Cat}(\mathbf{z}_n | \boldsymbol{\pi}) = \prod_{k=1}^K \pi_k^{z_{nk}} $$


 $\mathbf{x}_n$と$\mathbf{z}_n$の同時分布は、依存関係から次のようになる。

$$ p(\mathbf{x}_n, \mathbf{z}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) = p(\mathbf{x}_n | \mathbf{z}_n, \boldsymbol{\mu}) p(\mathbf{z}_n | \boldsymbol{\pi}) $$

 $\mathbf{x}_n$の周辺分布を考える。$\mathbf{x}_n$と$\mathbf{z}_n$の同時分布を$\mathbf{z}_n$に関して周辺化する。ただし、$\mathbf{z}_n$の全ての組み合わせ($(z_{n1}, z_{n2}, \cdots, z_{nK})$として$(1, 0, \cdots, 0)$から$(0, 0, \cdots, 1)$まで)の和を$\sum_{\mathbf{z}}$で表す。

$$ \begin{aligned} p(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) &= \sum_{\mathbf{z}_n} p(\mathbf{x}_n, \mathbf{z}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) \\ &= \sum_{\mathbf{z}_n} p(\mathbf{x}_n | \mathbf{z}_n, \boldsymbol{\mu}) p(\mathbf{z}_n | \boldsymbol{\pi}) \\ &= \sum_{\mathbf{z}_n} \prod_{k=1}^K p(\mathbf{x}_n | \boldsymbol{\mu}_k)^{z_{nk}} \prod_{k=1}^K \pi_k^{z_{nk}} \\ &= \sum_{\mathbf{z}_n} \prod_{k=1}^K \Bigl( \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \Bigr)^{z_{nk}} \end{aligned} $$

 $\sum_{z}$と$\prod_{k=1}^K$を展開すると

$$ \begin{aligned} p(\mathbf{x}) &= \sum_{\mathbf{z}_n} \Bigl\{ \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^{z_{n1}} \Bigl( \pi_2 p(\mathbf{x}_n | \boldsymbol{\mu}_2) \Bigr)^{z_{n2}} \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^{z_{nK}} \Bigr\} \\ &= \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^1 \Bigl( \pi_2 p(\mathbf{x}_n | \boldsymbol{\mu}_2) \Bigr)^0 \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^0 \\ &\qquad + \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^0 \Bigl( \pi_2 p(\mathbf{x}_n | \boldsymbol{\mu}_2) \Bigr)^1 \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^0 \\ &\qquad \vdots \\ &\qquad + \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^0 \Bigl( \pi_2 p(\mathbf{x}_n | \boldsymbol{\mu}_2) \Bigr)^0 \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^1 \\ &= \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) + \pi_2 p(\mathbf{x}_n | \boldsymbol{\mu}_2) + \cdots + \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \end{aligned} $$

となるので、1から$K$の和をまとめると

$$ p(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) = \sum_{k=1}^K \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \tag{9.47} $$

$\mathbf{x}_n$の周辺分布は混合ベルヌーイ分布の定義式になるのが分かる。

 以上で混合ベルヌーイ分布を確認できた。次からは、混合ベルヌーイ分布に対する最尤推定を考える。

・不完全データ対数尤度の確認

 $N$個の観測データ$\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_N\}$が独立して生成されると仮定すると、不完全データ尤度関数は次のように分解できる。

$$ \begin{aligned} p(\mathbf{X} | \boldsymbol{\mu}, \boldsymbol{\pi}) &= \prod_{n=1}^N p(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) \\ &= \prod_{n=1}^N \sum_{k=1}^K \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \\ &= \prod_{n=1}^N \sum_{k=1}^K \pi_k \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) \end{aligned} $$

 よって、不完全データ対数尤度関数は次の式になる。

$$ \begin{align} \ln p(\mathbf{X} | \boldsymbol{\mu}, \boldsymbol{\pi}) &= \sum_{n=1}^N \ln p(\mathbf{x}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) \\ &= \sum_{n=1}^N \ln \left\{ \sum_{k=1}^K \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \right\} \tag{9.51}\\ &= \sum_{n=1}^N \ln \left\{ \sum_{k=1}^K \pi_k \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) \right\} \end{align} $$

 不完全データ対数尤度関数は$\ln$の中に$\sum$を含むため解析的に解けない。そこで、完全データ対数尤度を扱う(潜在変数を導入する)。

・完全データ対数尤度の確認

 観測データ$\mathbf{X}$に加えて潜在変数$\mathbf{Z} = \{\mathbf{z}_1, \mathbf{z}_2, \cdots, \mathbf{z}_N\}$が与えられている(完全データ)とすると、完全データ対数尤度関数は次のように分解できる。

$$ \begin{aligned} p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\mu}, \boldsymbol{\pi}) &= \prod_{n=1}^N p(\mathbf{x}_n, \mathbf{z}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) \\ &= \prod_{n=1}^N p(\mathbf{x}_n | \mathbf{z}_n, \boldsymbol{\mu}) p(\mathbf{z}_n | \boldsymbol{\pi}) \\ &= \prod_{n=1}^N \prod_{k=1}^K p(\mathbf{x}_n | \boldsymbol{\mu})^{z_{nk}} \prod_{k=1}^K \pi_k^{z_{nk}} \\ &= \prod_{n=1}^N \prod_{k=1}^K \Bigl\{ \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}) \Bigr\}^{z_{nk}} \\ &= \prod_{n=1}^N \prod_{k=1}^K \left\{ \pi_k \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) \right\}^{z_{nk}} \\ &= \prod_{n=1}^N \prod_{k=1}^K \left\{ \pi_k \prod_{i=1}^D \mu_{ki}^{x_{ni}} \ln (1 - \mu_{ki})^{(1 - x_{ni})} \right\}^{z_{nk}} \end{aligned} $$

 よって、完全データ対数尤度関数は次の式になる。

$$ \begin{align} \ln p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\mu}, \boldsymbol{\pi}) &= \sum_{n=1}^N \ln p(\mathbf{x}_n, \mathbf{z}_n | \boldsymbol{\mu}, \boldsymbol{\pi}) \\ &= \sum_{n=1}^N \sum_{k=1}^K z_{nk} \Bigl\{ \ln \pi_k + \ln p(\mathbf{x}_n | \boldsymbol{\mu}) \Bigr\} \\ &= \sum_{n=1}^N \sum_{k=1}^K z_{nk} \left[ \ln \pi_k + \sum_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) \right] \\ &= \sum_{n=1}^N \sum_{k=1}^K z_{nk} \left[ \ln \pi_k + \sum_{i=1}^D \Bigl\{ x_{ni} \ln \mu_{ki} + (1 - x_{ni}) \ln (1 - \mu_{ki}) \Bigr\} \right] \tag{9.54} \end{align} $$

 しかし、実際には潜在変数$\mathbf{Z}$は観測できないので、$\mathbf{Z}$の事後分布による期待値を用いることにする。$\mathbf{Z}$の事後分布は、観測データ$\mathbf{X}$から求められる。

・対数尤度の期待値の導出

 潜在変数$\mathbf{Z}$の事後分布は、尤度を式(9.52)、事前分布を式(9.53)、周辺分布を(9.47)として、ベイズの定理を用いて求められる。

$$ \begin{aligned} p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\pi}) &= \frac{ p(\mathbf{X} | \mathbf{Z}, \boldsymbol{\mu}) p(\mathbf{Z} | \boldsymbol{\pi}) }{ \sum_{\mathbf{Z}} p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\mu}, \boldsymbol{\pi}) } \\ &= \prod_{n=1}^N \frac{ p(\mathbf{x}_n | \mathbf{z}_n, \boldsymbol{\mu}) p(\mathbf{z}_n | \boldsymbol{\pi}) }{ \sum_{\mathbf{z}_n} p(\mathbf{x}_n | \mathbf{z}_n, \boldsymbol{\mu}) p(\mathbf{z}_n | \boldsymbol{\pi}) } \\ &= \prod_{n=1}^N \frac{ \prod_{k=1}^K \Bigl\{ \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}) \Bigr\}^{z_{nk}} }{ \sum_{\mathbf{z}_n} \prod_{k=1}^K \Bigl\{ \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}) \Bigr\}^{z_{nk}} } \\ &= \prod_{n=1}^N \frac{ \prod_{k=1}^K \Bigl\{ \pi_k \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) \Bigr\}^{z_{nk}} }{ \sum_{\mathbf{z}_n} \prod_{k=1}^K \Bigl\{ \pi_k \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) \Bigr\}^{z_{nk}} } \end{aligned} $$

 この分布を用いて、完全データ対数尤度の期待値を求める。

 対数尤度の期待値は、離散確率分布の期待値の定義式(1.33)より、次の式で求められる。

$$ \begin{aligned} \sum_{\mathbf{Z}} p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\pi}) \ln p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\mu}, \boldsymbol{\pi}) &= \sum_{\mathbf{Z}} p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\mu}, \boldsymbol{\pi}) \sum_{n=1}^N \sum_{k=1}^K z_{nk} \left[ \ln \pi_k + \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \right] \\ \Leftrightarrow \mathbb{E}_{p(\mathbf{Z} | \mathbf{X})} [ \ln p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\mu}, \boldsymbol{\pi}) ] &= \mathbb{E}_{p(\mathbf{Z} | \mathbf{X})} \left[ \sum_{n=1}^N \sum_{k=1}^K z_{nk} \Bigl[ \ln \pi_k + \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \Bigr] \right] \\ &= \sum_{n=1}^N \sum_{k=1}^K \mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n)} [ z_{nk} ] \left[ \ln \pi_k + \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \right] \end{aligned} $$

 この式に含まれる$\mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n)} [z_{nk}]$を求める。

$$ \begin{aligned} \mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n)} [ z_{nk} ] &= \sum_{\mathbf{z}_n} z_{nk} p(\mathbf{z}_n | \mathbf{x}_n, \boldsymbol{\mu}, \boldsymbol{\pi}) \\ &= \frac{ \sum_{\mathbf{z}_n} z_{nk} \prod_{k'=1}^K \Bigl\{ \pi_{k'} p(\mathbf{x}_n | \boldsymbol{\mu}_{k'}) \Bigr\}^{z_{nk'}} }{ \sum_{\mathbf{z}_n} \prod_{k'=1}^K \Bigl\{ \pi_{k'} p(\mathbf{x}_n | \boldsymbol{\mu}_{k'}) \Bigr\}^{z_{nk'}} } \end{aligned} $$

 この式の分母は、先ほど確認した周辺分布(9.47)より、混合分布に置き換えられる。

$$ \begin{aligned} \mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n)} [ z_{nk} ] &= \frac{ \sum_{\mathbf{z}_n} z_{nk} \prod_{k'=1}^K \Bigl\{ \pi_{k'} p(\mathbf{x}_n | \boldsymbol{\mu}_{k'}) \Bigr\}^{z_{nk'}} }{ \sum_{k'=1}^K \pi_{k'} p(\mathbf{x}_n | \boldsymbol{\mu}_{k'}) } \end{aligned} $$

 また分子は、$\sum_{\mathbf{z}_n}$と$\prod_{k'=1}^K$を展開すると

$$ \begin{aligned} \sum_{\mathbf{z}_n} z_{nk} \prod_{k'=1}^K \Bigl\{ \pi_{k'} p(\mathbf{x}_n | \boldsymbol{\mu}_{k'}) \Bigr\}^{z_{nk'}} &= \sum_{\mathbf{z}_n} z_{nk} \Bigl\{ \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^{z_{n1}} \cdots \Bigl( \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \Bigr)^{z_{nk}} \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^{z_{nK}} \Bigr\} \\ &= 0 * \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^1 \cdots \Bigl( \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \Bigr)^0 \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^0 \\ &\qquad \vdots \\ &\qquad + 1 * \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^0 \cdots \Bigl( \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \Bigr)^1 \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^0 \\ &\qquad \vdots \\ &\qquad + 0 * \Bigl( \pi_1 p(\mathbf{x}_n | \boldsymbol{\mu}_1) \Bigr)^0 \cdots \Bigl( \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \Bigr)^0 \cdots \Bigl( \pi_K p(\mathbf{x}_n | \boldsymbol{\mu}_K) \Bigr)^1 \\ &= \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) \end{aligned} $$

となり(たぶん?)、$k' = k$の因子のみが残る。

$$ \begin{align} \mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n)} [ z_{nk} ] &= \frac{ \pi_k p(\mathbf{x}_n | \boldsymbol{\mu}_k) }{ \sum_{k'=1}^K \pi_{k'} p(\mathbf{x}_n | \boldsymbol{\mu}_{k'}) } \\ &= \frac{ \pi_k \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{ki}) }{ \sum_{k'=1}^K \pi_{k'} \prod_{i=1}^D \mathrm{Bern}(x_{ni} | \mu_{k'i}) } = \gamma(z_{nk}) \tag{9.56} \end{align} $$

 観測データによって学習した$z_{nk} = 1$となる期待値$\mathbb{E}_{p(\mathbf{z}_n | \mathbf{x}_n)} [z_{nk}]$を負担率と呼び、$\gamma(z_{nk})$とおく。

 $\gamma(z_{nk})$に置き換えた式を$L$とおく。

$$ \mathbb{E}_{p(\mathbf{Z} | \mathbf{X})} [ \ln p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\mu}, \boldsymbol{\pi}) ] = \sum_{n=1}^N \sum_{k=1}^K \gamma(z_{nk}) \left[ \ln \pi_k + \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \right] \equiv L \tag{9.55} $$

 $L$を最大化する各パラメータを求めていく。

・パラメータの最尤解の導出

 完全データ対数尤度の期待値$L$を最大化するパラメータ$\boldsymbol{\mu}$を求めていく。

 $L$を$\mu_{ki}$に関して微分する。

$$ \begin{aligned} \frac{\partial L}{\partial \mu_{ki}} &= \sum_{n=1}^N \left[ \frac{\partial}{\partial \mu_{ki}} \sum_{k=1}^K \gamma(z_{nk}) \ln \pi_k + \frac{\partial}{\partial \mu_{ki}} \sum_{k=1}^K \gamma(z_{nk}) \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \right] \\ &= \sum_{n=1}^N \gamma(z_{nk}) \frac{\partial}{\partial \mu_{ki}} \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \end{aligned} $$

 $D$個のベルヌーイ分布の微分は

$$ \begin{aligned} \frac{\partial}{\partial \mu_{ki}} \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) &= \frac{\partial}{\partial \mu_{ki}} \sum_{i=1}^D \Bigl\{ x_{ni} \ln \mu_{ki} + (1 - x_{ni}) \ln (1 - \mu_{ki}) \Bigr\} \\ &= \frac{\partial}{\partial \mu_{ki}} \sum_{i=1}^D x_{ni} \ln \mu_{ki} + \frac{\partial}{\partial \mu_{ki}} \sum_{i=1}^D (1 - x_{ni}) \ln (1 - \mu_{ki}) \\ &= x_{ni} \frac{1}{\mu_{ki}} + (1 - x_{ni}) \frac{1}{(1 - \mu_{ki})} (- 1) \end{aligned} $$

となる。自然対数の微分$(\ln x)' = \frac{1}{x}$と、$g(x) = 1 - \mu_{ki}$、$f(g(x)) = \ln (1 - \mu_{ki})$として合成関数の微分$\{f(g(x))\}' = f'(g(x)) g'(x)$を行った。
 上の式に代入すると、完全データ対数尤度の期待値の微分$\frac{\partial}{\partial \mu_{ki}}$が得られる。

$$ \begin{aligned} \frac{\partial}{\partial \mu_{ki}} &= \sum_{n=1}^N \gamma(z_{nk}) \left\{ x_{ni} \frac{1}{\mu_{ki}} - (1 - x_{ni}) \frac{1}{(1 - \mu_{ki})} \right\} \end{aligned} $$

 $\frac{\partial}{\partial \mu_{ki}}$を0とおき、$\mu_{ki}$に関して解く。

$$ \begin{aligned} \frac{\partial}{\partial \mu_{ki}} &= \sum_{n=1}^N \gamma(z_{nk}) \left\{ x_{ni} \frac{1}{\mu_{ki}} - (1 - x_{ni}) \frac{1}{(1 - \mu_{ki})} \right\} = 0 \\ &\Leftrightarrow \sum_{n=1}^N \gamma(z_{nk}) x_{ni} \frac{1}{\mu_{ki}} - \sum_{n=1}^N \gamma(z_{nk}) (1 - x_{ni}) \frac{1}{(1 - \mu_{ki})} = 0 \\ &\Leftrightarrow \sum_{n=1}^N \gamma(z_{nk}) x_{ni} (1 - \mu_{ki}) - \sum_{n=1}^N \gamma(z_{nk}) (1 - x_{ni}) \mu_{ki} = 0 \\ &\Leftrightarrow \mu_{ki} \sum_{n=1}^N \gamma(z_{nk}) = \sum_{n=1}^N \gamma(z_{nk}) x_{ni} \\ &\Leftrightarrow \mu_{ki} = \frac{1}{N_k} \bar{x}_{ni} \end{aligned} $$

 $\mu_{ki}$の最尤解が得られた。ただし

$$ \begin{align} N_k &= \sum_{n=1}^N \gamma(z_{nk}) \tag{9.57}\\ \bar{x}_{ki} &= \frac{1}{N_k} \sum_{n=1}^N \gamma(z_{nk}) x_{ni} \tag{9.58} \end{align} $$

とおく。
 他の$i = 1, 2, \cdots, D$の項も同様に求められるので

$$ \bar{\mathbf{x}}_k = (\bar{x}_{k1}, \bar{x}_{k2}, \cdots, \bar{x}_{kD})^{\top} $$

とおくと、クラスタ$k$のパラメータ$\boldsymbol{\mu}_k$の最尤解が得られる。

$$ \boldsymbol{\mu}_k = \bar{\mathbf{x}}_k \tag{9.59} $$

 他のクラスタについても同様に求められる。

・混合係数の最尤解の導出

 最後に、完全データ対数尤度の期待値$L$を最大化する混合係数を求めていく。

 $\boldsymbol{\pi}$には総和が1という制約条件(9.9)があるため、ラグランジュ未定乗数法を用いて$L$に制約条件を含めた式を立てて$F(\boldsymbol{\pi})$とおく。

$$ \begin{aligned} F(\boldsymbol{\pi}) &= L + \lambda \left( \sum_{k=1}^K \pi_k - 1 \right) \\ &= \sum_{n=1}^N \sum_{k=1}^K \gamma(z_{nk}) \left\{ \ln \pi_k + \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \right\} + \lambda \left( \sum_{k=1}^K \pi_k - 1 \right) \end{aligned} $$

 $F(\boldsymbol{\pi})$を$\pi_k$に関して微分する。

$$ \begin{aligned} \frac{\partial F(\boldsymbol{\pi})}{\partial \pi_k} &= \sum_{n=1}^N \left\{ \frac{\partial}{\partial \pi_k} \sum_{k=1}^K \gamma(z_{nk}) \ln \pi_k + \frac{\partial}{\partial \pi_k} \sum_{k=1}^K \sum_{i=1}^D \ln \mathrm{Bern}(x_{ni} | \mu_{ki}) \right\} + \frac{\partial}{\partial \pi_k} \lambda \left( \sum_{k=1}^K \pi_k - 1 \right) \\ &= \sum_{n=1}^N \gamma(z_{nk}) \frac{1}{\pi_k} + \lambda \end{aligned} $$

 $\frac{\partial F(\boldsymbol{\pi})}{\partial \pi_k}$は9.2.2項の式(9.21)と同じ式になる。
 したがって、$\frac{\partial F(\boldsymbol{\pi})}{\partial \pi_k}$を0とおき、$\pi_k$に関して解くと、9.2.2項「混合ガウス分布のEMアルゴリズム」と同じ手順で最尤解が得られる。

$$ \pi_k = \frac{1}{N} \sum_{n=1}^N \gamma(z_{nk}) \tag{9.60} $$

 他のクラスタについても同様に求められる。

 以上でパラメータ$\boldsymbol{\mu},\ \boldsymbol{\pi}$の最尤解を求められた。

参考文献

  • C.M.ビショップ著,元田 浩・他訳『パターン認識と機械学習 上下』,丸善出版,2012年.

おわりに

 ブログにまとめているのは文章の練習のためでもあるので、本によって解説のレベル感や文体を変えているのですが、この本はまだ定まっていない。本当はですます調で書きたい。

 実装編を書いている間に更新式の意味についての理解が深まった(けどこの数理編には書き足せていない)ので、冒頭にリンクを張った記事も参考にしてもらうとイメージしやすくなるかもしれません。

 2021年5月28日は、アンジュルムの佐々木莉佳子さんの20歳の誕生日です!おめでとうございます!

 先日公開された新曲です♪最初と最後に映っている方です。