EM 演算法

一個簡單的範例

假如我們觀察拋擲銅板的現象，得到觀察序列 x = {0, 1, 0, 0, 1, 1, 0, 0, 0, 1 } 這個現象，其中的 1 代表正面 (人頭)，0 代表反面 (字)，因此正面共出現 4 次，反面共出現 6 次。

因此，p’(1) = 0.4, p’(0) = 0.6。

那麼我們應該如何假設 p(0) 與 p(1) 的機率分布呢？

根據最大似然法則，我們應該去找出一個機率模型 p 可以最大化下列算式。

\arg\max_p \; p(x) = \arg\max \prod_i p(x_i)

我們可以計算看看下列兩個機率模型 p1, p2, p3 可能產生 x 的機率各為多少。

\begin{aligned} p1(1)=0.5 ; p1(0) = 0.5 \\ p2(1)=0.2 ; p2(0) = 0.8 \\ p3(1)=0.4 ; p3(0) = 0.6\end{aligned}

根據簡單的機率公式，我們可以算出下列結果。

\begin{aligned} p1(x) = \prod_i p1(x_i) = 0.5^4 * 0.5^6 = 0.00097656 \\ p2(x) = \prod_i p2(x_i) = 0.2^4 * 0.8^6 = 0.00041943 \\ p3(x) = \prod_i p3(x_i) = 0.4^4 * 0.6^6 = 0.00119439\end{aligned}

因此果然驗證了最可能的機率模型是 p3，也就是 p(1)=0.4, p(0)=0.6。

雖然我們找出了符合觀察現象 x 的最可能機率模型 (p3) ，但是對於投擲銅板這件事而言，p3 卻不是最適當的模型，因為最適當的模型是 p1 。

這個例子說明了一件事實，用最大似然法則所找出來的機率模型 p’ 未必是真正的機率源模型，只是根據觀察現象 x 所推導出來的最佳化機率模型而已。

但是，假如統計資料 x 序列的長度更長，那麼 x 的統計數據通常會更接近真實機率分布 X，因此最大似然法則所找出的機率模型 p’ 也就會更接近機率源模型 p，於是我們就可以認為 p’ 足以代表 p 了。

銅板問題的最大概似估計

EM 演算法是一種「最大概似估計」 (Maximum Likelihood Estimation, MLE)，要瞭解 EM 演算法之前，先讓我們瞭解何謂「最大概似估計」。

假如連續投擲一個銅板，結果有 H 次正面，T 次反面，那麼假設該銅板的正面機率為 $\theta$ ，那麼請問甚麼樣的 $\theta$ 會讓這個 (H, T) 結果的機率最大呢？

對於正面機率 $\theta$ 的銅板，我們可以用二項分布計算出現 #T 次正面 #H 次反面的機率為 $C(H+T, H) * \theta^H (1-\theta)^T$ 。

由於其中的 C(H+T, H) 與 $\theta$ 無關，因此我們只需要最大化後面那一項，也就是 $\arg\max_{\theta} \theta^H (1-\theta)^T$ 。

此時若我們先取 log，則最大值的 $\theta$ 點並不會改變，因此我們可以改為 $\arg\max_{\theta}\;H \log(\theta)+T\log(1-\theta)$

到底甚麼樣的 $\theta$ 值會讓上述算式最大呢？我們可以對上式取微分，尋找斜率為零的點。

$\frac{d}{d\theta} (H \log(\theta)+T\log(1-\theta)) = \frac{H}{\theta}+\frac{T}{1-\theta} = 0$ ; 連續可微函數最大值的斜率等於零。

=> $(1-\theta) H = T \theta$ ; 移項可得

=> $H=\theta (T+H)$ ; 將 $\theta$ 放到同一邊

=> $\theta = \frac{H}{T+H}$ ; 求得解答，最大化該式的 $\theta$ 為 $`\frac{H}{T+H}$$

因此，我們可以知道最大化 $\arg\max_{\theta}\;\theta^H (1-\theta)^T$ 這個算式的 $\theta$ 為 $\frac{H}{T+H}$ ，而這個 $\theta$ 也正是該問題的最大概似估計。

最大條件機率的分布

針對許多機率現象，我們只能觀察到某些面向的結果，但是無法觀察到全部的面向。這種情況就可以使用條件機率。

根據最大似然法則，假如已觀察到聯合機率分布 (X,Y)，其中 (x,y) 事件出現的機率為 p’(x,y) ，那麼根據最大似然法則，我們應當尋求盡可能滿足下列條件的算式。

\arg\max_h \; P(x,y|h)

然而，通常雙變數的聯合機率分布 p’(x,y) 會遭遇到『樣本稀疏性』的問題，因此若直接最大化上述算式，將會造成相當大的統計偏差。

為了解決『樣本稀疏性』的問題，我們應該採用較為可信的 p’(x) 作為 p(x) 的估計，p’(y) 作為 p(y) 的估計，而非直接採用 p’(x,y) 作為 p(x,y) 的估計值。

p'(x,y) = \frac{p'(x,y)}{p'(x) p'(y)} * p'(x) p'(y) \sim \frac{p(x,y)}{p'(x) p'(y)} * p'(x) p'(y)

於是我們可以最佳化下列算式

\arg\max_p \; \frac{p(x,y)}{p'(x) p'(y)} * p'(x) p'(y)

根據條件機率的定義，我們可以將 p’(x,y) 改寫如下。

p'(x,y) = p'(x)*p'(y|x)

如果我們用 p(y|x) 取代 p’(y|x)，那麼我們應該最大化下列算式。

於是我們可以最大化下列算式。

\arg\max \; { p'(x) * p(y|x) }

針對機率分布 p 而言，其機率為 p(x,y) 相當於下列算式。

arg\max \; p(X',Y') = \arg\max \; \prod_{x',y'} p(x',y') =  \arg\max \; \prod_{x',y'} p(x') p(y'|x')

根據微積分的原理，如果我們對上述算式進行微分的動作，那麼最佳解將會式微分式為 0 的 p 解 。