ＥＭアルゴリズムは、クラスター分析の一種で教師無し学習に分類されます。単純なk-平均法とは違い、どのクラスターに属しているか確率を考慮しています。
仮定
・（入力データ）データを確率変数$x$で表します。
・（クラスターの数）データをＫ個に分類。Ｋは定数で分析者が決めます。
・（隠れ変数）データがどのクラスターに属しているかを表した確率変数$z_i$（i=1~K）を導入します。属しているクラスターの番号のときだけ１で他は０をとるとします。さらに下記条件を満たすと仮定します。
条件１．$P(z_i=1)=\pi_i$,$P(z_i=0)=1-\pi_i$
条件２．$\sum_{i=1}^{K}\pi_i=1$
$z_i$たちを$K$までまとめた確率変数の組み合わせを$z=(z_1,・・・,z_K)$と書くことにすると、あるデータが2つ目のクラスターに属しているときは$z=(0,1,0,0,・・・,0)$となります。
・（クラスターごとの分布）属するクラスターが決まったときの、データの密度関数$f_i(x|z_i=1)$は平均$μ_i$、分散$\Sigma_i$の正規分布$N(x|μ_i,\Sigma_i)$に従うとします。ただしこの正規分布の次元はデータの特徴量の次元と等しいとします。
アルゴリズムの準備
・$\pi_i$と$z$の仮定より$z$の確率関数は、$h(z)=\prod_{i=1}^K\pi_i^{z_i}$となります。$z$の中である１つの$z_j$だけ１になるので、$\pi_j$だけ残りこの等号が成立します。
・クラスターごとの分布の仮定より、$f(x|z)=\prod_{i=1}^{K}N(x|\mu_i,\Sigma_i)^{z_i}$となります。
・$x$と$z$の同時分布は代入して$g(x,z)=h(z)f(x|z)=\prod_{i=1}^K\pi_i^{z_i}$ ×$\prod_{i=1}^{K}N(x|μ_i,\Sigma_i)^{z_i}$となります。
・上記の周辺化によりデータ$x$の分布は$f(x)=\sum_{zのﾊﾟﾀｰﾝ全て} g(x,z)$であり、$z$のパターンは$z=(1,0,・・・,0)～(0,0,・・・,0,1)$しかないので、最終的に$$f(x)=\sum_{i=1}^K\pi_{i}N(x|\mu_i,\Sigma_i)$$となる。これを混合ガウス分布といいます。これでデータの分布モデルができました。
・$x$の分布を設定できたので、最尤法により実際の入力データ$x_1,・・・,x_N$が独立で同分布の密度関数$f(x)$に従うとして対数尤度$L(\pi,\mu,Σ)=\ln\prod_{n=1}^Nf(x_n)$をとると
$L(\pi,\mu,Σ)=\sum_{n=1}^N\ln f(x_n)$
$=\sum_{n=1}^N\ln \sum_{i=1}^K \pi_{i} N(x_n,\mu_i,\Sigma_i)$.
これを最大にするパラメータ$\pi_i$ , $\mu_i$ , $Σ_i(i=1～K)$を求めたいわけです。
パラメータの最適化
対数尤度にした後は、それぞれのパラメータで偏微分した式を０とおいて、連立方程式を解くことになります。この記事で計算は省略しますが、実際にそうして計算すると下記のようになります。
$\mu_i=\frac{\sum_{n=1}^Nγ(z_{ni})x_n}{\sum_{n=1}^Nγ(z_{ni})}$
$\Sigma_i=\frac{1}{\sum_{n=1}^Nγ(z_{ni})}\sum_{n=1}^Nγ(z_{ni})(x_n-\mu_i)(x_n-\mu_i)^{T}$
$\pi_i=\frac{\sum_{n=1}^Nγ(z_{ni})}{N}$
※ただし、$γ(z_{ni})=\frac{\pi_iN(x_n|\mu_i,\Sigma_i)}{\sum_{i=1}^K\pi_iN(x_n|\mu_i,\Sigma_i)}$ とおいています。これを負担率といいます。
このようにパラメータを計算しましたが、$\pi_i$ , $\mu_i$ , $\Sigma_i(i=1～K)$の全てに$γ(z_{ni})$が入っており、$γ(z_{ni})$自身が$\pi_i$ , $\mu_i$ , $\Sigma_i(i=1～K)$によって決まるため、堂々巡りとなります。そこで、先にまずランダムでパラメータ$\pi_i$ , $\mu_i$ , $\Sigma_i(i=1～K)$から決めてやって、そこから負担率を算出して、その負担率から２回目のパラメータを求めるといった繰り返し処理によって最適なパラメータに近づくようにします。
ＥＭアルゴリズム
１．（初期化）パラメータ$\pi_i$ , $\mu_i$ , $\Sigma_i(i=1～K)$を初期化する。
－ループ開始－
２．（負担率の計算）パラメータ$\pi_i$ , $\mu_i$ , $\Sigma_i(i=1～K)$を用いて負担率を計算する。
３．上記のパラメータ最適化の計算から、２の負担率を使ってパラメータ$\pi_i$ , $\mu_i$ , $\Sigma_i(i=1～K)$を算出する。
４．対数尤度の式に３で求めたパラメータを代入する。前回の対数尤度と比較して変化無しまたは設定しておいた差よりも小さくなったところでアルゴリズムを終了する。そうでなければ２に戻る。
－ループ終了－
◆大事な事
・パラメータを更新するごとに最適解に近づくことが数学的に知られています。
・条件付き分布は正規分布でなくてもよく、もっと一般的な理論が展開されています。
参考
・はじめてのパターン認識（平井 有三）
・http://seetheworld1992.hatenablog.com/entry/2018/10/01/002958