PRMLの第9章に、EMアルゴリズムの適応対象として、混合ベルヌーイ分布(Bernoulli Mixture Models)を紹介しています。混合ベルヌーイ分布は、EMアルゴリズムを勉強する際に非常に有用な応用事例なので、今回はMINISTの画像クラスタリングを混合ベルヌーイ分布で行い、その学習を通してEMアルゴリズムについて理解を深めていきましょう。
MNISTデータセットは、28×28サイズの白黒画像で、合計784ピクセルの画像になります。
今回は、混合ベルヌーイ分布を用いて、これらのMNISTデータセットのクラスタリングをやっていきます。
- 混合ベルヌーイ分布を用いて、MNISTデータセットのクラスタリングを行う
- EMアルゴリズムでパラメータの最適化を行う
混合ベルヌーイ分布によるMNISTクラスタリングの全体像
今回、この記事でやりたいことである、混合ベルヌーイ分布を用いたクラスタリングの全体像のアルゴリズムをまとめます。

MNISTは上記のような、28×28= 784ピクセルからなる、手書き文字のデータセットです。
そのため1枚の手書き文字画像の次元数を$D$とすると、$D = 784$となります。今回、混合ベルヌーイ分布を用いて、これらの手書き文字のクラスタリングをやっていきます。
MNISTで用意されている手書き文字は、0~9の10種類あるので、今回クラスタリングによって分けるクラス数を$K$とすると $K=10$となり、それぞれの要素は $K = {1, 2, \dots , 10 }$となります。
MNISTデータセットの用意
まず、MNISTデータセットを準備するところから始めましょう。手元にMNISTのデータセットがすぐに用意できる人はそれを用いて構いません。
今回はsklearnを用いて、MNISTのデータセットを用意します。
また今回は混合ベルヌーイ分布を用いて手書き数字のクラスタリングに取り組みますが、ベルヌーイ分布は0と1の2値しか表現できないため、0~255の256値になっているMNISTのデータセットを前処理する必要性があります。
scikit-learnを用いてMNISTデータセットをダウンロード
詳しいMNISTの説明やデータセットのダウンロード方法は、こちらを参考にしてください。

MNISTデータセットのダウンロード
混合ベルヌーイ分布の定式化
まず、1枚のMNIST画像をベルヌーイ分布を用いて表現することを考えていきます。
今、MNISTは28×28の$D = 784$ ピクセルのデータです。
ここで、各ピクセルは0,1 の二値画像を想定しているので、各ピクセルは下記のベルヌーイ分布で表現することができます。
2値画像のあるピクセル $x_i $( $i = \{1, 2, \dots, D\} $) は、パラメータ$\mu_i$を持つ、下記のベルヌーイ分布で表現できる
\begin{equation} p(x_i) = Bern(x_i | \mu_i) = \mu_i^{x_i}(1-\mu_i)^{1-x_i} \end{equation}
784ピクセルある画像のうち、1ピクセルはこのように表現することができました。では次に、これを全ピクセル$D=784$次元に拡張してみましょう。各ピクセルが独立のベルヌーイ分布にしたがっているとすると、$D$次元の2値画像の同時確率分布は下記のように表現することができます。
ベクトル$\bm{x}$を、$ \bm{x} = \{ x_1, x_2, \dots, x_D \}$とすると、D次元の2値画像の確率分布は、D個のベルヌーイ分布を用いて次のように表現することができる。
\begin{equation} p(\bm{x} | \bm{\mu}) = \prod_{i=1}^{D} \mu_i^{x_i}(1-\mu_i)^{1-x_i} \end{equation}
単純に、(1)式をD個分掛け合わせたのが(2)式です。
ここで、確率変数$\bm{x}$については、その期待値はベルヌーイ分布の期待値から、次のようになります。
\mathbb{E}[\bm{x}] = \bm{\mu}
ここで、$\bm{\mu}$は、$\bm{\mu} = \{ \mu_1, \mu_2, \dots, \mu_D \}$である。
続いて、(2)式で表すことができた、D次元のピクセル画像の確率$p(\bm{x})$ですが、これらの有限混合分布を考えてみることにします。
混合ベルヌーイ分布を用いてK=10のクラスタリング
混合ガウス分布の例でもあるように、今回は$K=10$種類の手書き文字をクラスタリングしたいと考えています。まず、クラス$k$に属する画像の分布を表現してみます。
その前に、(2)式をそのまま用いて、クラスkに所属する画像の確率分布を表現してみます。
\begin{equation} p(\bm{x} | \bm{\mu_k}) = \prod_{i=1}^{D} \mu_{ki}^{x_i}(1-\mu_{ki})^{1-x_i} \end{equation}
ここで、$\bm{mu_k}$は、$D$次元のパラメータベクトルになります。
通常の正規分布では、パラメータは2個ですが、今回の混合ベルヌーイ分布では、ピクセルの数$D$だけパラメータがあることになります(!)
では、これらの分布からなる有限混合分布はこのようになります。
混合比率を$\bm{\pi} = \{ \pi_1, \pi_2, \dots, \pi_K \}$とすると、次のように表現することができます。
\begin{equation} p(\bm{x} | \bm{\mu}, \bm{\pi}) = \sum_{k=1}^{K} \pi_k p(\bm{x} | \bm{\mu_k}) \end{equation}
ここで、$\bm{\mu}$は下記のベクトルである。
\bm{\mu} = \{ \bm{\mu_1}, \bm{\mu_2}, \dots, \bm{\mu_k}\}
ただし、$\bm{\pi} = \{ \pi_1, \pi_2, \dots, \pi_K \}$の各要素$\pi_k$は下記の条件を満たす。
\int_{k=1}^{K} \pi_k = 1
ここの、有限混合分布の表現までは混合ガウス分布(GMM)の例とほとんど一緒なので、比較してみると、混合分布が何をやろうとしているかについて、ざっくりとでも理解することができるだろう。
混合ベルヌーイ分布の尤度関数を考える
いよいよ、EMアルゴリズムで学習をする前段階に入っていきます。
今、データ$N$個のデータ $\bm{X} = \{x_1, x_2, \dots, x_N\}$が与えれているとします。MNISTのデータセットの場合、訓練データが70,000個あるため、$N=70,000$です。
(4)式にこれらのデータを全て入れた尤度関数$\mathcal{L}$を考えると、
\begin{equation} \begin{split} \mathcal{L} &= \prod_{n=1}^{N} p(\bm{x_n} | \bm{\mu}, \bm{\pi}) \\ &= \prod_{n=1}^{N} \biggr \{ \sum_{k=1}^{K} \pi_{k} p(\bm{x} | \bm{\mu_k}) \biggr\} \end{split} \end{equation}
となるので、例のごとく対数尤度関数$ln\mathcal{L}$を考えるため、(5)式の尤度関数にlogをとると
\begin{equation} \begin{split} ln\mathcal{L} &= ln \biggr \{ \prod_{n=1}^{N} p(\bm{x_n} | \bm{\mu}, \bm{\pi}) \biggr \} \\ &= \sum_{n=1}^{N} ln \biggr \{ \sum_{k=1}^{K} \pi_{k} p(\bm{x} | \bm{\mu_k}) \biggr\} \end{split} \end{equation}
ここまでで、モデル化し、対数尤度関数まで得ることができました。
機械学習におけるモデルの学習法の1つとして、データが得られているときは(6)のような対数尤度関数(6を最大化するようなパラメータを見るける、最尤推定法がありますが、(6)式を眺めると、logの中身に$\sum$にが入った形式になっております。
これはいわゆる、ログサムなどと呼ばれており、残念ながら解析的に解を得ることができないことが知られています。
実際、これらの尤度関数を最大化するようなパラメータを見るけるためには、EMアルゴリズムを利用します。
よくテキストでは、混合ガウスモデル(GMM)でEMアルゴリズムの説明などで登場するようなEMアルゴリズムですが、今回のような混合ベルヌーイモデルにおいてもEMアルゴリズムがこの尤度関数の最大化(つまり最尤推定)を行うのに利用することができます。
EMアルゴリズムで推論
潜在変数$\bm{z}$の導入
まず、EMアルゴリズムを利用するために、潜在変数$\bm{z}$を考えます。
この潜在変数$\bm{z}$は、1-of-K表現になっており、$\bm{z} = \{ z_1, z_2, …, z_k\}$ かつ、$z \in {0, 1}$かつ、$ \sum_{n=1}{k} z_n = 1$となっています。
隠れ変数zは下記のような分布になっている。
\bm{z} = \{ z_1, z_2, ..., z_k\}
z_i \in {0, 1} \\ \sum_{k=1}^{K} z_k = 1
この隠れ変数を導入して、(4)式の分布を考えてみましょう。(4)式では、$\bm{\pi}$を用いましたが、隠れ変数を導入すると、$\bm{\pi}$の代わりに$\bm{z}$で、表現することになります。
(4)式の確率分布は、隠れ変数$\bm{z} $を用いると、このように書くことができます。
\begin{equation} \begin{split} p(\bm{x} | \bm{z}, \bm{\mu} ) &= \prod_{k=1}^{K} p(\bm{x} | \bm{\mu_k})^{z_k} \end{split} \end{equation}
ここで、潜在変数$\bm{z}$の事前分布は下記のように表現することができます。
\begin{equation} p(\bm{z}| \bm{\pi}) = \prod_{k=1}^{K} \pi_{k}^{z_k} \end{equation}
ここで、EMアルゴリズムを考えるために、完全データ$p(\bm{x}, \bm{z}|\bm{\pi}, \bm{\mu})$の分布を考えていきます。
(7)式と(8)式を用いると、条件付き確率の連鎖律の公式から、下記のような式変形をすることで、完全データの分布を得ることができます。
\begin{equation} \begin{split} p(\bm{x}, \bm{z}|\bm{\pi}, \bm{\mu}) &= p(\bm{x} | \bm{z}, \bm{\mu} ) p(\bm{z}| \bm{\pi}) \\ &= \prod_{k=1}^{K} p(\bm{x} | \bm{\mu_k})^{z_k} \prod_{k=1}^{K} \pi_{k}^{z_k} \\ &= \prod_{k=1}^{K} \biggr ( \pi_k p(\bm{x} | \bm{\mu_k})^{z_k} \biggr )^{z_k} \\ &= \prod_{k=1}^{K} \biggr ( \pi_k \biggr ( \prod_{i=1}^{D} \mu_{ki}^{x_i}(1-\mu_{ki})^{1-x_i} \biggr ) \biggr )^{z_k} \\ \end{split} \end{equation}
と表現することができます。
最後の式変形は、(3)の式を代入して行いました。
続いて、EMアルゴリズムを行うために、完全データ対数尤度関数を考えます。
完全データ対数尤度関数は、下記のような尤度関数である。
\begin{equation} \begin{split} ln p(\bm{X}, \bm{Z}|\bm{\pi}, \bm{\mu}) &= ln \biggr ( \prod_{n=1}^{N} p(\bm{x_n}, \bm{z_n}|\bm{\pi}, \bm{\mu}) \biggr ) \end{split} \end{equation}
ここで、$\bm{X}$は各データセット、$\bm{Z}$は完全データのベクトルです。
(10)式を変形すると、下記のような完全データ対数尤度を得ることができます。
\begin{equation} \begin{split} ln p(\bm{X}, \bm{Z}|\bm{\pi}, \bm{\mu}) &= \sum_{n=1}^{N} ln p(\bm{x_n}, \bm{z_n} | \bm{\pi},\bm{\mu}) \\ &= \sum_{n=1}^{N} ln \biggr \{ \prod_{k=1}^{K} \biggr ( \pi_k \biggr ( \prod_{i=1}^{D} \mu_{ki}^{x_{ni}}(1-\mu_{ki})^{1-x_{ni}} \biggr ) \biggr )^{z_{nk}} \biggr \} \\ &= \sum_{n=1}^{N} \sum_{k=1}^{K} z_{nk} \biggr \{ ln \pi_k + \sum_{i=1}^{D} \biggr ( x_{ni} ln \mu_{ki}+ (1-x_{ni}) ln (1-\mu_{ki}) \biggr ) \biggr \} \\ \end{split} \end{equation}
これで、完全データ対数尤度関数を得ることができた。
EMアルゴリズムでの学習
(11)式で、完全データ対数尤度関数を得ることができたので、ここからはEMアルゴリズムを用いてパラメータの最適化をやっていきます。
Eステップにおける学習
Eステップにおける学習では、潜在変数$\bm{z}$の事後分布$\bm{Z}$を求めることになります。
EMアルゴリズムのEステップでは、潜在変数の事後分布を求めることになる。