混合ガウス分布の最尤推定を解く手法のひとつとしてEMアルゴリズムを導入する．

混合ガウス分布

KK 個のクラスがあり，それぞれを zkz_k と表記する．そしてその実現値を1-of-K符号化で z\boldsymbol{z} と表す(どれか一つだけが1，残りが0)．各クラスの周辺分布は p(z)=πkp(\boldsymbol{z}) = \pi_k であるとする．そしてクラスが既知である場合の条件付き確率は

p(x∣zk=1)=N(x∣μk,Σk)\begin{align*} p (\boldsymbol{x} \mid z_{k}=1) = \mathcal{N} (\boldsymbol{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}) \end{align*}

であるとすると， x\boldsymbol{x} の周辺分布は

p(x)=∑k=1KπkN(x∣μk,Σk)\begin{align*} p(\boldsymbol{x})=\sum_{k=1}^{K} \pi_{k} \mathcal{N} (\boldsymbol{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}) \end{align*}

となる．逆に x\boldsymbol{x} が与えられた場合の z\boldsymbol{z} の条件付き確率は

γ(zk)=p(zk=1∣z)=p(zk=1)p(x∣zk=1)∑j=1Kp(zj=1)p(x∣zj=1)=πkN(x∣μk,Σk)∑j=1KπjN(x∣μj,Σj)\begin{align*} \gamma (z_k) = p(z_k = 1 \mid \boldsymbol{z}) &= \dfrac{p(z_k = 1)p(\boldsymbol{x} \mid z_k = 1)}{\sum_{j=1}^{K}p(z_j = 1)p(\boldsymbol{x} \mid z_j = 1)} \\ &= \dfrac{\pi_{k} \mathcal{N} (\boldsymbol{x} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}(\boldsymbol{x} \mid \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j})} \end{align*}

で与えられる．これは x\boldsymbol{x} がクラスkを説明する度合いであると言え，負荷率と呼ばれる．

最尤推定の特異性

GMMから得られた観測値の集合 {x1,…,xN}\{ \boldsymbol{x}_1, \ldots, \boldsymbol{x}_N \} について X\mathbf{X} を各行が xiT\boldsymbol{x}_{i}^{\text{T}} である行列を X\mathbf{X} とする．すると尤度関数の対数は

ln⁡p(X∣π,μ,Σ)=∑n=1Nln⁡{∑k=1KπkN(xn∣μk,Σk)}\begin{align*} \ln p(\mathbf{X} \mid \boldsymbol{\pi}, \boldsymbol{\mu}, \mathbf{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N} (\boldsymbol{x}_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k})\right\} \end{align*}

で与えられる．ここではシンプルに Σk=σk2I\mathbf{\Sigma}_k = \sigma_{k}^{2} \mathbf{I} とし，訓練データのうちの一つ xn\boldsymbol{x}_n がある要素の平均 μj\boldsymbol{\mu}_j と等しかったとする．このとき指数関数部が1になるので上の式の ln⁡\ln の和の中に

N(xn∣xn,σj2I)=1(2π)D/21σjD\begin{align*} \mathcal{N}(\mathbf{x}_{n} \mid \mathbf{x}_{n}, \sigma_{j}^{2} \mathbf{I})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{\sigma_{j}^{D}} \end{align*}

の項が現れる．補足 で見るように指数がかかった項は σ→0\sigma \rightarrow 0 で0に収束するが，この項は無限大に発散してしまう．よってこの対数尤度関数は無限大に発散し，尤度関数も発散する．

補足

1次元のガウス分布に従う x∼N(μ,σ2)x \sim \mathcal{N}(\mu, \sigma^2) から得られたデータ [x1,…,xN][x_1, \ldots, x_N] のうち xix_i が μ\mu に一致しているとする．これに対して尤度関数を求めると

12πσ∏n≠iN12πσe−(xn−μ)22σ2\begin{align*} \frac{1}{\sqrt{2\pi} \sigma} \prod_{n \neq i}^N \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{(x_n-\mu)^2}{2\sigma^2}} \end{align*}

となる．ここで σ→0\sigma \rightarrow 0 とすると指数部の収束が 1/σ1 / \sigma の発散より強いため0に収束する．

GMMのEMアルゴリズム

対数尤度を μk\boldsymbol{\mu}_k について微分することで以下を得る．

γnk=πkN(xn∣μk,Σk)∑jπkN(xn∣μj,Σj)Nk=∑n=1Nγnk\begin{align*} \gamma_{nk} &= \dfrac{\pi_k \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k)}{\sum_j \pi_k \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_j, \mathbf{\Sigma}_j)} \\ N_k &= \sum_{n=1}^{N} \gamma_{nk} \end{align*}

さらにこの値を使って

μk‾=1Nk∑n=1NγnkxnΣk‾=1Nk∑n=1Nγnk(xn−μk‾)(xn−μk‾)Tπk‾=NkN\begin{align*} \overline{\boldsymbol{\mu}_k} &= \dfrac{1}{N_k} \sum_{n=1}^{N} \gamma_{nk} \boldsymbol{x}_n \\ \overline{\mathbf{\Sigma}_{k}} &= \dfrac{1}{N_{k}} \sum_{n=1}^{N} \gamma_{nk} (\boldsymbol{x}_{n} - \overline{\boldsymbol{\mu}_{k}}) (\boldsymbol{x}_{n} - \overline{\boldsymbol{\mu}_{k}})^{\text{T}} \\ \overline{\pi_k} &= \dfrac{N_k}{N} \end{align*}

を得る．このように得られた値で (μk,Σk,πk)(\boldsymbol{\mu}_k, \mathbf{\Sigma}_k, \pi_k) を更新する処理を繰り返す．

EMアルゴリズムの一般形

観測値 X={x1,…,xN}\mathbf{X} = \{ \boldsymbol{x}_1, \ldots, \boldsymbol{x}_N \} が潜在変数(例えばどのGMMにおいてどのクラスに属するか，など)を含む確率分布 p(X,Z∣θ)p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta} ) に従うとする．これは人間が決めるものである．そして p(X∣θ)p( \mathbf{X} \mid \boldsymbol{\theta} ) を最大化するような θ\boldsymbol{\theta} を求める(ここで X\mathbf{X} は観測値なので定数，よって θ\boldsymbol{\theta} の関数になる)．

EMアルゴリズムでは潜在変数の事後分布 p(Z∣X,θ)p( \mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta} ) を用いる．これも以下のようにして

p(Z∣X,θ)=p(X,Z∣θ)∫p(X,Z∣θ)dZ\begin{align*} p( \mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta} ) = \frac{p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{\int p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}) d \mathbf{Z}} \end{align*}

前提で仮定した同時分布 p(X,Z∣θ)p( \mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta} ) のみから求められる．初期推定値 θ0\boldsymbol{\theta}_0 を用いて以下の処理を繰り返す．

  1. p(Z∣X,θi)p( \mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}_i ) を計算する
  2. Q(θ)Q(\boldsymbol{\theta}) を ∫p(Z∣X,θ)ln⁡p(X,Z∣θi)dZ\int p (\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta} ) \ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}_i) d\mathbf{Z} とする
  3. この最大値を与える θ\boldsymbol{\theta} を θi+1\boldsymbol{\theta}_{i+1} とする

GMMのEMアルゴリズム では負荷率 γnk\gamma_{nk} を求めた部分が(1)の事後分布を求める部分に相当する．この枠組みに当てはめると p(Z∣xn,θold)p( \mathbf{Z} \mid \boldsymbol{x}_n, \boldsymbol{\theta}_{\mathrm{old}} ) は

p(Z∣xn,θ)=∏k=1K(πkN(xn∣μk,Σk)∑jπjN(xn∣μj,Σj))znk\begin{align*} p( \mathbf{Z} \mid \boldsymbol{x}_n, \boldsymbol{\theta} ) = \prod_{k=1}^{K} \left( \dfrac{\pi_k \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_j, \mathbf{\Sigma}_j)} \right)^{z_{nk}} \end{align*}

で与えられる．同様に p(xn,Z∣θ)p(\boldsymbol{x}_n, \mathbf{Z} \mid \boldsymbol{\theta}) は

p(xn,Z)=∏k=1KπkznkN(xn∣μk,Σk)znk\begin{align*} p( \boldsymbol{x}_n, \mathbf{Z}) = \prod_{k=1}^{K} \pi_{k}^{z_{nk}}\mathcal{N}( \boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k)^{z_{nk}} \end{align*}

で与えられるから対数は

ln⁡p(xn,Z)=∑k=1Kznkln⁡πk+znkln⁡N(xn∣μ,Σk)\begin{align*} \ln p( \boldsymbol{x}_n, \mathbf{Z}) = \sum_{k=1}^{K} {z_{nk}} \ln \pi_{k} + {z_{nk}} \ln \mathcal{N}( \boldsymbol{x}_n \mid \boldsymbol{\mu}, \mathbf{\Sigma}_k) \end{align*}

となる．ここで Z\mathbf{Z} は1-of-K符号化により {(1,0,…,0),(0,1,0,…,0),…,(0,0,…,1)}\{ (1, 0, \ldots ,0), (0, 1, 0, \ldots, 0), \ldots, (0, 0, \ldots, 1) \} を動くとして

∫p(Z∣xn,θold)ln⁡p(xn,Z∣θ)dZ=∑k=1Kγnk{ln⁡πk+ln⁡N(xn∣μk,Σk)}\begin{align*} \int p (\mathbf{Z} \mid \boldsymbol{x}_n, \boldsymbol{\theta}_{\mathrm{old}} ) \ln p(\boldsymbol{x}_n, \mathbf{Z} \mid \boldsymbol{\theta}) d\mathbf{Z} = \sum_{k=1}^{K} \gamma_{nk} \{ \ln \pi_k + \ln \mathcal{N}(\boldsymbol{x}_n \mid \boldsymbol{\mu}_k, \mathbf{\Sigma}_k) \} \end{align*}

となる．これを xn\boldsymbol{x}_n につき和をとったのが(2)の Q(θ)Q(\boldsymbol{\theta}) に相当する． (μk,Σk,πk)(\boldsymbol{\mu}_k, \mathbf{\Sigma}_k, \pi_k) について微分した連立方程式を0にする値を求めるとGMMのEMアルゴリズム で求めた式が得られる．

そもそも潜在変数とは

個人的には説明変数として用いている x\boldsymbol{x} 以外の切り捨てられた変数全てのことなのかなと思った．実際説明変数を列挙しようとすれば，いくらでもきりがなく出てくると思われる．あるいはハイパーパラメータと言ってもよいかもしれない．

例

ベイズ線形回帰

TODO

RVM

TODO

EMアルゴリズムが対数尤度を極大化する理由

一般的な議論を行う．以下の尤度関数

p(X∣θ)=∑Zp(X,Z∣θ)\begin{align*} p(\mathbf{X} \mid \boldsymbol{\theta})=\sum_{\mathbf{Z}} p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}) \end{align*}

を θ\boldsymbol{\theta} について最大化したい．

ln⁡p(X∣θ)=L(q,θ)+KL(q∣p)\begin{align*} \ln p(\mathbf{X} \mid \boldsymbol{\theta})=\mathcal{L}(q, \boldsymbol{\theta})+\mathrm{KL}(q | p) \end{align*}

ここで KL(q∣p)\mathrm{KL}(q | p) はKLダイバージェンスと呼ばれるが常に0以上の値をとる．以下の図は対数尤度を分解した図である．

対数尤度の分解

Eステップでは q(Z)q(\mathbf{Z}) について L(q,θ)\mathcal{L}(q, \boldsymbol{\theta}) を最大化する．すると下の図のように KL(q∣p)\mathrm{KL}(q | p) が0になる．

q(Z)に関する最大化(Eステップ)

次のMステップでは L(q,θ)\mathcal{L}(q, \boldsymbol{\theta}) を θ\boldsymbol{\theta} について最大化する．このとき KL(q∣p)\mathrm{KL}(q | p) は必ず0以上であるから，図のように ln⁡(X∣θ)\ln ( \mathbf{X} \mid \boldsymbol{\theta}) も点線から増加する．

Mステップ