EM算法(Expectation Maximization Algorithm)

1. 前言

??這是本人寫的第一篇博客(2013年4月5日發在cnblogs上,現在遷移過來)，是學習李航老師的《統計學習方法》書以及斯坦福機器學習課Andrew Ng的EM算法課后，對EM算法學習的介紹性筆記，如有寫得不恰當或錯誤的地方，請指出，并多多包涵，謝謝。另外本人數學功底不是很好，有些數學公式我會說明的仔細點的，如果數學基礎好，可直接略過。

2.基礎數學知識

??在正式介紹EM算法之前，先介紹推導EM算法用到的數學基礎知識，包括凸函數，Jensen不等式。

2.1.凸函數

??對于凸函數，凹函數，如果大家學過高等數學，都應該知道，需要注意的是國內教材如同濟大學的《高等數學》的這兩個概念跟國外剛好相反，為了能更好的區別，本文章把凹凸函數稱之為上凸函數，下凸函數，具體定義如下：

上凸函數：函數 $f(x)$ 滿足對定義域上任意兩個數 $a$ , $b$ 都有 $f[(a+b)/2] ≥ [f(a)+f(b)]/2$
下凸函數：函數 $f(x)$ 滿足對定義域上任意兩個數 $a$ , $b$ 都有 $f[(a+b)/2] ≤ [f(a)+f(b)]/2$

??更直觀的可以看圖2.1和2.2：


圖2.1. 上凸函數	圖2.2. 下凸函數

??可以清楚地看到圖2.1上凸函數中，$f[(a+b)/2] ≥ [f(a)+f(b)]/2$，而且不難發現，如果f(x)是上凸函數，那么 $-f(x)$ 是下凸函數。
??當 $a≠b$ 時，$f[(a+b)/2] > [f(a)+f(b)]/2$ 成立，那么稱 $f(x)$ 為嚴格的上凸函數，等號成立的條件當且僅當 $a=b$,下凸函數與其類似。

2.2.Jensen不等式

??有了上述凸函數的定義后，我們就能很清楚的Jensen不等式的含義，它的定義如下：

如果f是上凸函數，$X$ 是隨機變量，那么 $f(E[X]) ≥ E[f(X)]$
特別地，如果f是嚴格上凸函數，那么 $E[f(X)] = f(E[X])$ 當且僅當 $p(X=E[X])=1$，也就是說 $X$ 是常量。

??那么很明顯 $\log x$ 函數是上凸函數，可以利用這個性質。
??有了上述的數學基礎知識后，我們就可以看具體的EM算法了。

3.EM算法所解決問題的例子

??在推導EM算法之前，先引用《統計學習方法》中EM算法的例子：

例1. (三硬幣模型)假設有3枚硬幣,分別記作 $A,B,C$ 。這些硬幣正面出現的概率分別為 $π$，$p$ 和 $q$。投幣實驗如下，先投 $A$，如果 $A$ 是正面，即 $A=1$，那么選擇投 $B$；$A=0$，投 $C$。最后，如果 $B$ 或者 $C$ 是正面，那么 $y=1$；是反面，那么 $y=0$；獨立重復 $n$ 次試驗 $(n=10)$，觀測結果如下： $1,1,0,1,0,0,1,0,1,1$ 假設只能觀測到投擲硬幣的結果，不能觀測投擲硬幣的過程。問如何估計三硬幣正面出現的概率，即 $\pi$，$p$ 和 $q$ 的值。

解：設隨機變量 $y$ 是觀測變量，則投擲一次的概率模型為:$$P(y|\theta)=\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y}$$有 $n$ 次觀測數據 $Y$，那么觀測數據 $Y$ 的似然函數為:$$P(Y|\theta) = \prod_n^{j=1}[\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]$$那么利用最大似然估計求解模型解，即
\begin{align}
\widehat{\theta}& = \arg \max_{\theta} \log P(Y|\theta)\label{ex:loglikelihood1} \\
& = \arg \max_{\theta} \sum_{j=1}^{10} \log P(y^j|\theta)\label{ex:loglikelihood2} \\
& = \arg \max_{\theta} \sum_{j=1}^{10} \log [\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]\label{ex:loglikelihood3}
\end{align} 這里將概率模型公式和似然函數代入 $\eqref{ex:loglikelihood1}$ 式中，可以很輕松地推出 $\eqref{ex:loglikelihood1} \Rightarrow \eqref{ex:loglikelihood2} \Rightarrow \eqref{ex:loglikelihood3}$，然后選取 $\theta(\pi,p,q)$，使得 $\eqref{ex:loglikelihood3}$ 式值最大，即最大似然。然后，我們會發現因為 $\eqref{ex:loglikelihood3}$ 中右邊多項式 $+$ 符號的存在，使得 $\eqref{ex:loglikelihood3}$ 直接求偏導等于 $\theta$ 或者用梯度下降法都很難求得 $\theta$ 值。
這部分的難點是因為 $\eqref{ex:loglikelihood3}$ 多項式中 $+$ 符號的存在，而這是因為這個三硬幣模型中，我們無法得知最后得結果是硬幣 $B$ 還是硬幣 $C$ 拋出的這個隱藏參數。那么我們把這個latent 隨機變量加入到 log-likelihood 函數中，得
\begin{align}
l(\theta)& = \sum_{j=1}^{10}log\: \sum_{i=1}^{2}P(y_{j},z_{i}\mid \theta )\label{ex:latentlog1} \\
& = \sum_{j=1}^{10}log\: \sum_{i=1}^{2}Q_j(z_{i})\frac{P(y_{j},z_{i}\mid \theta )}{Q_j(z_{i})}\label{ex:latentlog2} \\
& \geq \sum_{j=1}^{10} \sum_{i=1}^{2}Q_j(z_{i})log\:\frac{P(y_{j},z_{i}\mid \theta )}{Q_j(z_{i})}\label{ex:latentlog3}
\end{align}
略看一下，好像很復雜，其實很簡單，首先是公式 $\eqref{ex:latentlog1}$，這里將 $z_i$ 做為隱藏變量，當 $z_1$ 為結果由硬幣 $B$ 拋出，$z_2$ 為結果由硬幣C拋出，不難發現:$$\sum_{i=1}^{2}P(y_{j},z_{i}\mid \theta )=P(y_{j}\mid \theta )
=\pi p^{y_{j}}(1-p)^{1-y_{j}}+\pi q^{y_{j}}(1-q)^{1-y_{j}}$$ ??接下來公式說明 $\eqref{ex:latentlog1} \Rightarrow \eqref{ex:latentlog2}$ (其中 $\eqref{ex:latentlog2}$ 中 $Q(z)$ 表示的是關于 $z$ 的某種分布，$\sum_iQ_j(z_i)=1$，很直接，在 $P$ 的分子分母同乘以 $Q_(z_i)$。最后是 $\eqref{ex:latentlog2} \Rightarrow \eqref{ex:latentlog3}$，到了這里終于用到了第二節介紹的Jensen不等式，數學好的人可以很快發現，$\sum_{i=1}^2Q_j(z_i)\frac{P(y_j,z_i|\theta)}{Q_j(z_i)}$ 就是 $[\frac{P(y_j,z_i|\theta)}{Q_j(z_i)}]$ 的期望值（如果不懂，可google期望公式并理解），且log是上凸函數，所以就可以利用Jensen不等式得出這個結論。因為我們要讓log似然函數 $l(\theta)$最大，那么這里就要使等號成立。根據Jensen不等式可得，要使等號成立，則要使 $\frac{P(y_j,z_i|\theta)}{Q_j(z_i)} =c$ 成立。
??再因為$\sum_iQ_j(z_i)=1$，所以得$\sum_iP(y_j,z_i|\theta)=c$，$c$ 為常數，那么(這里感謝網友@無影隨想指出錯誤)$$Q(z_{i})=P(y_{j},z_{i}\mid \theta )/\sum_{i}P(y_{j},z_{i}\mid \theta )
=P(y_{j},z_{i})/P(y_{j}\mid \theta ) =P(z_{i}\mid y_{j},\theta)$$這里可以發現$$Q_j(z_{1}) =\frac{\pi p^{y_{j}}(1-p)^{1-y_{j}}}{\pi p^{y_{j}}(1-p)^{1-y_{j}}+(1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}}\\
Q_j(z_{2} ) =\frac{(1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}}{\pi p^{y_{j}}(1-p)^{1-y_{j}}+(1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}}$$ ??OK,到這里，可以發現公式 $\eqref{ex:latentlog3}$ 中右邊多項式已經不含有“+”符號了，如果知道 $Q(z)$ 的所有值，那么可以容易地進行最大似然估計計算，但是 $Q$ 的計算需要知道 $\theta$ 的值。這樣的話，我們是不是可以先對θ進行人為的初始化 $\theta_0$，然后計算出 $Q$ 的所有值 $Q_1$ (在 $\theta_0$ 固定的情況下，可在 $Q_1$ 取到公式 $\eqref{ex:latentlog3}$ 的極大值)，然后在對公式 $\eqref{ex:latentlog3}$ 最大似然估計，得出新的 $\theta_1$ 值（在固定Q1的情況下，取到公式 $\eqref{ex:latentlog3}$ 的極大值），這樣又可以計算新的 $Q$ 值 $Q_1$,然后依次迭代下去。答案當然是可以。因為 $Q_1$ 是在 $\theta_0$ 的情況下產生的，可以調節公式 $\eqref{ex:latentlog3}$ 中 $\theta$ 值，使公式 $\eqref{ex:latentlog3}$ 的值再次變大，而 $\theta$ 值變了之后有需要調節 $Q$ 使 $\eqref{ex:latentlog3}$ 等號成立，結果又變大，直到收斂（單調有界必收斂），如果到現在還不是很清楚，具體清晰更廣義的證明可以見下部分EM算法說明。
??另外對公式 $\eqref{ex:latentlog3}$ 進行求偏導等于 $\theta$，求最大值，大家可以自己練習試試，應該很簡單的，這里不做過多陳述。
??在《統計學習方法》書中，進行兩組具體值的計算

$\pi_0=0.5,\ p_0=0.5,\ q_0=0.5$，迭代結果為 $\pi=0.5,\ p=0.6,\ q=0.5$
$\pi_0=0.4,\ p_0=0.6,\ q_0=0.7$，迭代結果為 $\pi=0.4064,\ p=0.5368,\ q=0.6432$
兩組值的最后結果不相同，這說明EM算法對初始值敏感，選擇不同的初值可能會有不同的結果，只能保證參數估計收斂到穩定點。因此實際應用中常用的辦法就是選取多組初始值進行迭代計算，然后取結果最好的值。

??在進行下部分內容之前，還需說明下一個東西。在上面的舉例說明后，其實可以發現上述的解決方法跟一個簡單的聚類方法很像，沒錯，它就是K-means聚類。K-means算法先假定k個中心，然后進行最短距離聚類，之后根據聚類結果重新計算各個聚類的中心點，一次迭代，是不是很像，而且K-means也是初始值敏感，因此其實K-means算法也包含了EM算法思想，只是這邊EM算法中用P概率計算，而K-means直接用最短距離計算。所以EM算法可以用于無監督學習。在下一篇文章，我準備寫下典型的用EM算法的例子，高斯混合模型(GMM,Gaussian Mixture Model)。

4.EM算法

4.1.模型說明

??考慮一個參數估計問題，現有 ${y_1,y_2,…,y_n}$ 共 $n$ 個訓練樣本，需有多個參數 $\pi$ 去擬合數據，那么這個 $\log$ 似然函數是：$$l(\theta) = \sum_{j=1}^{n} \log P(y_j|\theta)$$??可能因為 $\theta$ 中多個參數的某種關系（如上述例子中以及高斯混合模型中的3類參數），導致上面的 $\log$ 似然函數無法直接或者用梯度下降法求出最大值時的 $\theta$ 值，那么這時我們需要加入一個隱藏變量 $z$，以達到簡化 $l(\theta)$,迭代求解 $l(\theta)$ 極大似然估計的目的。

4.2.EM算法推導

??這小節會對EM算法進行具體推導，許多跟上面例子的解法推導是相同的，如果已經懂了，可以加速閱讀。首先跟“三硬幣模型”一樣，加入隱變量 $z$ 后，假設 $Q(z)$ 是關于隱變量 $z$ 的某種分布，那么有如下公式：
\begin{align}
l(\theta)& = \sum_{j=1}^{n}log\: \sum_{i=1}P(y_{j},z_{i}\mid \theta )\label{infer1}\\
& = \sum_{j=1}^{n}log\: \sum_{i=1}Q(z_{i})\frac{P(y_{j},z_{i}\mid \theta )}{Q(z_{i})}\label{infer2}\\
& \geq \sum_{j=1}^{n} \sum_{i=1}Q(z_{i})log\:\frac{P(y_{j},z_{i}\mid \theta )}{Q(z_{i})}\label{infer3}
\end{align}??公式 $\eqref{infer1}$ 是加入隱變量，$\eqref{infer1} \Rightarrow \eqref{infer2}$ 是在基礎上分子分母同乘以，$\eqref{infer2} \Rightarrow \eqref{infer3}$ 用到Jensen不等式（跟“三硬幣模型”一樣），等號成立的條件是，$c$ 是常數。再因為，則有如下 $Q$ 的推導：
\begin{equation*}\sum_{i}P(y_{j},z_{i}\mid \theta )/c=1\\
\Rightarrow \sum_{i}P(y_{j},z_{i}\mid \theta )=c\\
\qquad \qquad \qquad \qquad \Rightarrow Q_{j}(z_{i})=P(y_{j},z_{i}\mid \theta )/\sum_{i}P(y_{j},z_{i}\mid \theta )\\
\qquad \qquad \qquad \qquad \qquad =P(y_{j},z_{i}\mid \theta )/P(y_{j}\mid \theta )\\
\qquad \qquad \qquad =P(z_{i}\mid y_{j},\theta )
\end{equation*}??再一次重復說明，要使 $\eqref{infer3}$ 等式成立，則 $Q_j(z_i)$ 為 $y_j,\ z$的后驗概率。算出 $Q_j(z_i)$ 后對 $\eqref{infer3}$ 就可以進行求偏導，以剃度下降法求得 $\theta$ 值，那么又可以計算新 $Q_j(z_i)$ 的值，依次迭代，EM算法就實現了。

EM 算法(1)
選取初始值 $\theta_0$ 初始化 $\theta$，$t=0$
Repeat {
??E步：
\begin{equation*}
\begin{split}
Q_{j}^{t}(z_{i})& =P(y_{j},z_{i}\mid \theta^{t} )/\sum_{i}P(y_{j},z_{i}\mid \theta^{t} )\\
& =P(y_{j},z_{i}\mid \theta^{t} )/P(y_{j}\mid \theta^{t} )\\
& =P(z_{i}\mid y_{j},\theta^{t} )
\end{split}
\end{equation*}??M步：
\begin{equation*}
\begin{split}
\theta^{t+1}& =arg\: max_{\theta }\: \sum_{j=1}^{n} \sum_{i}Q_{j}^{t}(z_{i})log\:\frac{P(y_{j},z_{i}\mid \theta )}{Q_{j}^{t}(z_{i})}\\
t& =t+1
\end{split}
\end{equation*}}直到收斂

4.3.EM算法收斂性證明

??當 $\theta$ 取到 $\theta_t$ 值時，求得$$\theta^{t+1} =arg\: max_{\theta }\: \sum_{j=1}^{n} \sum_{i}Q_{j}^{t}(z_{i})log\:\frac{P(y_{j},z_{i}\mid \theta )}{Q_{j}^{t}(z_{i})}$$那么可得如下不等式：\begin{align}
l(\theta^{t+1} )& = \sum_{j=1}^{n}log\: \sum_{i}Q_{j}^{t}(z_{i})\frac{P(y_{j},z_{i}\mid \theta^{t+1} )}{Q_{j}^{t}(z_{i})}\label{orderProof1}\\
& \geq \sum_{j=1}^{n}\sum_{i}Q_{j}^{t}(z_{i})log\: \frac{P(y_{j},z_{i}\mid \theta^{t+1} )}{Q_{j}^{t}(z_{i})}\label{orderProof2}
\\
& \geq \sum_{j=1}^{n}\sum_{i}Q_{j}^{t}(z_{i})log\: \frac{P(y_{j},z_{i}\mid \theta^{t} )}{Q_{j}^{t}(z_{i})}\label{orderProof3}
\end{align}??$\eqref{orderProof1} \Rightarrow \eqref{orderProof2}$是因為Jensen不等式，因為等號成立的條件是 $\theta$ 為 $\theta_t$ 的時候得到的，而現在中的 $\theta$ 值為 $\theta_{t+1}$，所以等號不一定成立，除非 $\theta_{t+1} = \theta_t$，
??$\eqref{orderProof2} \Rightarrow \eqref{orderProof3}$ 是因為 $\theta_{t+1}$ 已經使得 $\sum_{j=1}^{n}\sum_{i}Q_{j}^{t}(z_{i})log\: \frac{P(y_{j},z_{i}\mid \theta^{t} )}{Q_{j}^{t}(z_{i})}$ 取得最大值，那必然不會小于 $\eqref{orderProof3}$ 式。
??所以 $l(\theta)$ 在迭代下是單調遞增的，且很容易看出 $l(\theta)$ 是有上界的 (單調有界收斂) ，則EM算法收斂性得證。

4.4. EM算法E步說明

??上述EM算法描述，主要是參考Andrew NG教授的講義，如果看過李航老師的《統計方法學》，會發現里面的證明以及描述表明上有些許不同，Andrew NG教授的講義的說明（如上述）將隱藏變量的作用更好的體現出來，更直觀，證明也更簡單，而《統計方法學》中則將迭代之間θ的變化羅列的更為明確，也更加準確的描述了EM算法字面上的意思：每次迭代包含兩步：E步，求期望；M步，求極大化。下面列出《統計方法學》書中的EM算法，與上述略有不同：

EM算法 (1)：
選取初始值θ0初始化θ，t=0
Repeat {
??E步：
\begin{equation}
\begin{split}
\label{Estep}
H(\theta ,\theta ^{t})& =E_{z}[logP(Y,Z\mid \theta )\mid Y,\theta^{t}]\\
& =\sum_{z}P(Z\mid Y,\theta ^{t})logP(Y,Z\mid \theta )
\end{split}
\end{equation}??M步：$$\theta^{t+1} = arg\: max_{\theta } \: H(\theta ,\theta^{t})$$}直到收斂

??$\eqref{Estep}$ 式中，$Y={y_1,y_2,…,y_m},\ Z={z_1,z_2,…,z_m}$,不難看出將 $\eqref{infer3}$ 式中兩個 $\sum$ 對換，就可以得出 $\eqref{Estep}$ 式，而 $\eqref{Estep}$ 式即是關于分布 $z$ 的一個期望值，而需要求這個期望公式，那么要求出所有的EM算法 (1) 中E步的值，所以兩個表明看起來不同的EM算法描述其實是一樣的。

5.小結

??EM算法的基本思路就已經理清，它計算是含有隱含變量的概率模型參數估計，能使用在一些無監督的聚類方法上。在EM算法總結提出以前就有該算法思想的方法提出，例如HMM中用的Baum-Welch算法就是。
??在EM算法的推導過程中，最精妙的一點就是 $\eqref{orderProof1}$ 式中的分子分母同乘以隱變量的一個分布，而套上了Jensen不等式，是EM算法順利的形成。