EM Algorithm 詳盡介紹: 利用簡單例子輕鬆讀懂 EM 的原理及概念

前言

在 Estimation – ML & MAP 這篇文章中我們介紹了 maximum likelihood (ML) 的概念，我們會先假設數據服從某一機率分布，在根據實驗數據以及 maximum log-likelihood 的方式尋找最佳的參數，但是有時候我們假設的機率分布可能會存在某個我們無法取得的潛在特徵，導致單純利用 ML 無法估計出未知參數，而 EM algorithm 就是幫助我們在遇到這種情況時依舊有辦法估計出未知參數，在本篇文章我們將會延用 Estimation – ML & MAP 投擲硬幣的例子，還沒有看過那篇的讀者可以點開稍微看一下實驗內容以及一些參數的定義。

ML Estimation Review

延續投擲硬幣的例子，假設我們今天從市面上蒐集 100 枚硬幣，這 100 枚硬幣由兩間不同的工廠 $A$、$B$ 製造，兩間工廠做出來的硬幣有不同的公差導致擲出正面的機率不同，幸運的是我們可以得知每一枚硬幣是由哪一間工廠製造的，令 $z_i = \{A=1, B=0\}$ 代表第 $i$ 枚硬幣出自哪一間工廠，$x_i = \{H = 1, T = 0\}$ 代表第 $i$ 枚硬幣投擲出正面或反面，而$$
\begin{aligned}
&p(z_i) = \theta^{z_i} (1-\theta)^{1-z_i} \\
&p(x_i | z_i = 1) = \alpha^{x_i} (1-\alpha)^{1-x_i} \\
&p(x_i | z_i = 0) = \beta^{x_i} (1-\beta)^{1-x_i} \\
&p(x_i | z_i) = \bigg[\alpha^{x_i} (1-\alpha)^{1-x_i} \bigg]^{z_i} \bigg[ \beta^{x_i} (1-\beta)^{1-x_i} \bigg]^{1-z_i}
\end{aligned}
$$ 分別代表：硬幣是由 $A$、$B$ 工廠製造的機率、由 $A$ 工廠製造的硬幣擲出正面反面的機率、由 $B$ 工廠製造的硬幣擲出正面反面的機率。需要注意的是，在 Estimation – ML & MAP 文章中我們假設的機率都是一個條件機率基於未知的參數，舉例來說我們應該把硬幣是由 $A$、$B$ 工廠製造的機率寫成 $p(z_i | \theta)$，但為了書寫方便我們一律將未知的參數 $\theta, \alpha, \beta$ 省略在條件機率裡面。

假設我們投擲完 100 枚硬幣後，發現由 $A$ 工廠製造的硬幣有 60 枚，並投出了 36 個正面、24 個反面，由 $B$ 工廠製造的應必有 40 枚，並投出了 16 個正面、24 個反面，我們可以利用 log maximum likelyhood 估測 3 個未知參數：$$
\begin{aligned}
\ln p(\mathbf{x}, \mathbf{z}) &= \ln p(\mathbf{z}) p(\mathbf{x} |\mathbf{z}) \\
&= \sum_{i=1}^{100} \ln p(z_i) p(x_i |z_i) \\
&= \sum_{i=1}^{100} \ln \theta^{z_i} (1-\theta)^{1-z_i} \bigg[\alpha^{x_i} (1-\alpha)^{1-x_i} \bigg]^{z_i} \bigg[ \beta^{x_i} (1-\beta)^{1-x_i} \bigg]^{1-z_i} \\
&= 60 \ln \theta + 40 \ln (1-\theta) + 36 \ln \alpha + 24 \ln (1- \alpha) + 16 \ln \beta + 24 \ln (1-\beta)
\end{aligned}
$$ 其中最後一個等式跳過了一些步驟，讀者可以利用 100 枚硬幣投擲結果自行帶入數學式中，即可得出相同結果，最後我們分別對 $\theta, \alpha, \beta$ 求導，可以估計出未知參數：$$
\begin{aligned}
&\frac{\partial \ln p(\mathbf{x}, \mathbf{z})}{\partial \theta} = \frac{60}{\theta} + \frac{40}{1-\theta} = 0 \rightarrow \theta = 0.6 \\
&\frac{\partial \ln p(\mathbf{x}, \mathbf{z})}{\partial \alpha} = \frac{36}{\alpha} + \frac{24}{1-\alpha} = 0 \rightarrow \alpha = 0.6 \\
&\frac{\partial \ln p(\mathbf{x}, \mathbf{z})}{\partial \beta} = \frac{16}{\beta} + \frac{24}{1-\beta} = 0 \rightarrow \beta = 0.4 \\
\end{aligned}
$$ 代表市面上流通的硬幣有 $60\%$ 的機率是由 $A$ 工廠產生，$40\%$ 由 $B$ 工廠產生，且 $A$ 工廠產生的硬幣有 $60\%$ 的機率會投擲出正面，$B$ 工廠則是有 $40\%$ 的機率投出正面。

Introduction – EM Algorithm

但是現實生活中不是時常這麼幸運的，我們可能必須在缺少某些資訊的情況下做出估計，舉例來說，我們可能不知道 100 枚硬幣中有哪幾枚是 $A$ 工廠製造，哪幾枚是 $B$ 工廠製造，我們只知道投擲 100 枚硬幣投出 52 次正面，48 次反面，此時 $z_i$ 是未知的、被隱藏的，但是這個隨機變數卻可能會影響到我們最終估測，這種變數我們稱之為潛在變數 latent variable。Latent variable 的出現使得我們無法單純利用 ML、MAP 估計未知參數，因為我們無法計算出 $p(\mathbf{x}, \mathbf{z})$，此時 EM algorithm 就派上用場了，它可以幫助我們在缺少某些資訊的情況下，依然可以利用 ML 估計未知參數。

當我們定義的模型或機率分布存在 latent variable 時，EM algorithm 可以幫助我們找出 maximum likelihood 的解。

如上述，likelihood $p(\mathbf{x}, \mathbf{z})$ 我們無法求得，那我們就束手無策了嗎？其實 likelihood 無法求得，我們可以求出 likelihood 的期望值，根據每個硬幣的正反結果 $x_i$，我們可以推估是由哪間工廠製作的機率，即 $p(z_i | x_i)$，再利用這個條件機率可以推出 likelihood 的期望值： $$
\mathbb{E}[\ln p(\mathbf{x}, \mathbf{z}) ] = \sum_{i=1}^{100} \bigg\{p(z_i=1 | x_i) \ln p(x_i, z_i=1) + p(z_i=0 | x_i) \ln p(x_i, z_i=0) \bigg\}
$$ 因為我們並不知道實際的 $z_i$ 是多少，導致 likelihood 未知，只能透過事後機率 $p(z_i | x_i)$ 計算 likelihood 的期望值，我們可以將 $p(z_i | x_i)$ 看做是一個權重：

• 當 $z_i = 1$ 時給定 $\ln p(x_i, z_i = 1)$ 權重 $p(z_i = 1| x_i)$
• 當 $z_i = 0$ 時給定 $\ln p(x_i, z_i = 0)$ 權重 $p(z_i = 0| x_i)$

藉此求出 likelihood 可能為多少，即其期望值 (平均值)。事後機率 $p(z_i | x_i)$ 的數學式可以由 Bayes’ theorem求出：$$
\begin{aligned}
p(z_i | x_i) = \frac{p(z_i) p(x_i | z_i)}{\sum_{z_i = 0}^1 p(z_i) p(x_i | z_i)} &= \frac{\theta^{z_i}(1-\theta)^{1-z_i} \bigg[\alpha^{x_i} (1-\alpha)^{1-x_i} \bigg]^{z_i} \bigg[ \beta^{x_i} (1-\beta)^{1-x_i} \bigg]^{1-z_i}} {\sum_{z_i = 0}^1 \theta^{z_i}(1-\theta)^{1-z_i} \bigg[\alpha^{x_i} (1-\alpha)^{1-x_i} \bigg]^{z_i} \bigg[ \beta^{x_i} (1-\beta)^{1-x_i} \bigg]^{1-z_i}} \\
\Rightarrow \hspace{0.5cm} p(z_i = 1 | x_i = 1) & = M = \frac{\theta\alpha}{\theta \alpha + (1-\theta) \beta } \\
\Rightarrow \hspace{0.5cm} p(z_i = 0 | x_i = 1) &= 1-M= \frac{(1-\theta) \beta}{\theta \alpha + (1-\theta) \beta } \\
\Rightarrow \hspace{0.5cm} p(z_i = 1 | x_i = 0) & = N= \frac{\theta (1-\alpha)}{\theta (1-\alpha) + (1-\theta) (1-\beta) } \\
\Rightarrow \hspace{0.5cm} p(z_i = 0 | x_i = 0) & = 1 – N = \frac{(1-\theta)(1-\beta)}{\theta (1-\alpha) + (1-\theta) (1-\beta) } \\
\end{aligned}
$$ 為了接下來的數學式簡潔，我們利用 $M$、$N$ 各自代表不同的機率值，將在接下來的推導使用。另外，由先前的推導 (ML Estimation Review) 我們亦可以求出 $p(x_i, z_i = 1)$ 以及 $p(x_i, z_i =0)$：$$
\begin{aligned}
p(x_i, z_i=1) &= \theta \alpha^{x_i} (1-\alpha)^{1-x_i} \\
p(x_i, z_i=0) &= (1-\theta) \beta^{x_i} (1-\beta)^{1-x_i} \\
\end{aligned}
$$

E step

至此我們可以開始 EM algorithm 了。第一步我們必須先有一個起始的參數值，好讓我們可以計算出事後機率，假設我們讓 $\theta^{\text{old}} = 0.5$，$\alpha^{\text{old}} = 0.5$，$\beta^{\text{old}} = 0.5$，我們就可以求得：$$
\begin{aligned}
p(z_i = 1 | x_i = 1) &= M = 0.5 \\
p(z_i = 0 | x_i = 1) &= 1-M = 0.5 \\
p(z_i = 1 | x_i = 0) &= N = 0.5 \\
p(z_i = 0 | x_i = 0) &= 1-N = 0.5
\end{aligned}
$$ 注意這裡的參數有一個上標 $\text{old}$，我們用來表示事後機率 $p(z_i | x_i)$ 裡的參數，而 $p(x_i, z_i)$ 的參數則沒有上標 $\text{old}$ 藉此區隔開兩個機率的參數。照理來說，$p(x_i, z_i)$ 和 $p(z_i | x_i)$ 裡的參數應該是同一個，但是 EM algorithm 為了估計出 likelihood 期望值，必須取得事後機率的實際值，才會將兩者的參數分開，具體原因以及實際操作會在下一步驟解釋。透過事後機率 $p(z_i | x_i)$ 初始的參數值，我們可以求出 $\mathbb{E}[\ln p(\mathbf{x}, \mathbf{z})]$，因此計算出期望值的這個步驟，我們稱為 expectation。$$
\begin{aligned}
\mathbb{E}[\ln p(\mathbf{x}, \mathbf{z})] &= \sum_{i=1}^{100} \bigg\{ p(z_i = 1|x_i) \ln \theta \alpha^{x_i} (1-\alpha)^{1-x_i} + p(z_i = 0|x_i) \ln (1-\theta) \beta^{x_i} (1-\beta)^{1-x_i} \bigg\} \\
&= 52M \ln \theta \alpha + 52 (1-M) \ln (1-\theta)\beta +\\
&\hspace{1cm} 48N \ln \theta( 1 – \alpha) + 48(1-N) \ln (1-\theta)(1-\beta)
\end{aligned}
$$
如上述，我們將參數起始值代入事後機率，$\theta^{\text{old}} = 0.5$，$\alpha^{\text{old}} = 0.5$，$\beta^{\text{old}} = 0.5$，藉此得到權重大小 $M$、$N$，但 $p(x_i, z_i)$ 的參數 $\theta$，$\alpha$，$\beta$ 則保留起來，並將於下一步利用偏微分對其做出估計。

M step

取得 likelihood 的期望值後，我們開始尋找此 likelihood 期望值最大值的參數點為何，作為估計結果，方法如同 ML、MAP，我們可以利用參數偏微分為 $0$ 的位置，找出最大值的點：
$$
\begin{aligned}
&\frac{\partial \mathbb{E}[\ln p(\mathbf{x}, \mathbf{z})]}{\partial \theta} = 0 \\
\Rightarrow \hspace{1cm} &\sum_{i=0}^{100} \bigg\{ \frac{p(z_i = 1 | x_i)}{\theta} – \frac{p(z_i = 0 | x_i)}{1-\theta} \bigg\} = 0 \\
\Rightarrow \hspace{1cm} &\frac{52M}{\theta} – \frac{52(1-M)}{(1-\theta)} + \frac{48N}{\theta} – \frac{48(1-N)}{(1-\theta)} = 0 \\
\Rightarrow \hspace{1cm} &\theta = \frac{52M + 48N}{100} = 0.5
\end{aligned}
$$
$$
\begin{aligned}
&\frac{\partial \mathbb{E}[\ln p(\mathbf{x}, \mathbf{z})]}{\partial \alpha} = 0 \\
\Rightarrow \hspace{1cm} &\sum_{i=0}^{100} \bigg\{ \frac{x_ip(z_i = 1 | x_i)}{\alpha} – \frac{(1-x_i)p(z_i = 0 | x_i)}{1-\alpha} \bigg\} = 0 \\
\Rightarrow \hspace{1cm} &\frac{52M}{\alpha} – \frac{48N}{1-\alpha} = 0 \\
\Rightarrow \hspace{1cm} &\alpha = \frac{52M}{52M + 48N} = 0.5 \\
\end{aligned}
$$
$$
\begin{aligned}
&\frac{\partial \mathbb{E}[\ln p(\mathbf{x}, \mathbf{z})]}{\partial \beta} = 0 \\
\Rightarrow \hspace{1cm} &\sum_{i=0}^{100} \bigg\{ \frac{x_ip(z_i = 1 | x_i)}{\beta} – \frac{(1-x_i)p(z_i = 0 | x_i)}{1-\beta} \bigg\} = 0 \\
\Rightarrow \hspace{1cm} &\frac{52(1-M)}{\beta} – \frac{48(1-N)}{1-\beta} = 0 \\
\Rightarrow \hspace{1cm} &\beta = \frac{52-52M}{100 – 52M – 48N} = 0.52 \\
\end{aligned}
$$

若讀者不清楚如何從第二式推到第三式，可以自行將 52 次正面 $x_i = 1$ 以及 48 次反面 $x_i=0$ 代入式子，搭配先前定義的 $M$、$1-M$、$N$、$1-N$ 等機率，即可推出相同的結果。此步驟涉及求出參數值使得期望值有最大值，因此我們將此步驟稱為 maximization。

EM algorithm

至此，從上述兩步驟我們可以一窺 EM algorithm 的概念，就是先給定一組優良的參數，算出事後機率 $p(z_i | x_i)$，在利用此事後機率計算出 likelihood 的期望值 $\mathbb{E}[\ln p(\mathbf{x}, \mathbf{z})]$，並對此期望值最偏微分取得最佳解的位置，接下來只要再次重複相同步驟直到收斂：$$
\{\theta^{0}, \alpha^{0}, \beta^{0}\} \rightarrow \text{E step: } \{M^{0}, N^{0}\} \rightarrow \text{M step: } \{\theta^{1}, \alpha^{1}, \beta^{1}\} \rightarrow \text{E step: } \{M^{1}, N^{1}\} \rightarrow …
$$

Discussion

透過以上述起始值以及 EM algorithm，迭代 10 次後我們可以求出收斂的 $\theta$，$\alpha$，$\beta$：$$
\begin{aligned}
\theta = 0.499468 \\
\alpha = 0.500000 \\
\beta = 0.539918
\end{aligned}
$$ 那如果我們從 $\theta = 0.7$，$\alpha = 0.5$，$\beta = 0.5$ 開始，迭代 10 次後我們會得到：$$
\begin{aligned}
\theta = 0.698567 \\
\alpha = 0.500000 \\
\beta = 0.564479
\end{aligned}
$$ 最後我們再試一組 $\theta = 0.6$，$\alpha = 0.59$，$\beta = 0.39$，迭代 10 次後得到：$$
\begin{aligned}
\theta = 0.604384 \\
\alpha = 0.590000 \\
\beta = 0.412943
\end{aligned}
$$ 從不同起始值出發我們可以得到不同的結果，但是這些參數值都有可能讓我們有相同的實驗結果，即投擲出 52 次正面、48 次反面，因此我們可以想見，我們擁有無限多種起始值，且這些起使值可能導致無限多種解，會有這種情況出現的原因是，我們缺少了 latent variable $z_i$ 的資訊，導致我們無法從這麼多解中取得一個最佳解。所以我們可以得到一個結論，雖然 EM algorithm 可以利用 likelihood 的期望值幫助我們估計出某個解，但是這個解會非常依賴參數的起始值，因此為了確保能夠找出不要過於偏離現實情況的解，我們必須對當前問題有基礎的了解，幫助我們判斷合理的起始值。

雖然可以利用 EM algorithm 找出具有 latent variable 之模型的解，但是該解非常依賴起始值，因此使用者最好對手邊問題有基礎了解，才能定出合理的初始值。

其他文章

機器學習系列文章總整理
Maximum Likelihood & Maximum a Posteriori－基礎估計模型的詳細介紹
Gaussian Mixture Model (GMM) －混合高斯模型的介紹以及完整數學式推導
[MobileNet 系列] MobileNetV1－為終端而生
[MobileNet 系列] MobileNetV2－小改動大進化，終端的唯一選擇