0x00 引言

EM算法是什么？什么是E（Epectation）？什么是M（Maximization）？什么又是公式里面出現的Q函數？這些公式都是怎么推導的？Nature拋硬幣的那個圖怎么就看不懂嘞？為什么看了那么多文章之后還是不懂？公式的符號怎么又不一樣呢？誰誰還說有九層塔？Emmm…interesting

下面，讓我們走進科學。

0x10 直觀理解

現在有一個隨機變量數據集$\mathbf{D}$，假設我們知道這個隨機變量$\mathbf{X}$服從某種分布(一般是高斯正規分布），我們的目的是想知道這個分布的參數$\theta$，可是隨機變量$\mathbf{X}$里面包含不可知的參數（也就是隱變量$\mathbf{z}$）的時候，EM算法是在一邊猜隱變量$\mathbf{z}$一邊更新$\theta$：

先蒙一下目標參數${\theta }^{old}$

E步：利用測算的目標參數${\theta }^{old}$和數據集$\mathbf{D}$猜隱變量$\mathbf{z}$的分布

M步：利用上一步猜出來的隱變量$\mathbf{z}$反思更新目標參數${\theta }^{new}$

重復上面兩步直到目標參數$\theta$收縮為止

再壓縮成人話的話，數據集$\mathbf{D}$是以$\mathbf{z}$為比例而測出來的$\mathbf{X}$，先用$\theta$蒙一個$\mathbf{z}$，然后用$\mathbf{z}$再算一個$\theta$，如此反復。

0x20 拋硬幣的例子

來源：Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm? Nature Biotechnology, 26(8), 897–899. https://doi.org/10.1038/nbt1406

上面這篇著名的EM入門論文里面有一張很好的圖例，利用拋硬幣來說明EM，可是對于某些初學者來講缺乏解讀可能還是有點難理解思路。

下面嘗試拆解一下分步驟解讀

0x21 問題定義

已知：

• 手上有兩種不同的硬幣，分別稱為A和B

實驗：

• 隨機拋硬幣十次為一組，記錄正面朝上（H）和反面朝上（T）的數據
• 換硬幣重復試驗

問題：

• 分別求這兩個硬幣正面朝上的概率${\theta }_A$和${\theta }_B$

?

0x22 完全信息 vs 包含隱函數的不完全信息

上圖的實驗過程中如果記錄了當時拋的是A或者B哪種硬幣，統計推斷的時候知道了每一組是屬于哪一種硬幣的情況下那當然很好算，這種情況叫完全信息。

假如實驗中根本不知道拋的時候究竟是哪一種硬幣，或者就不告訴你的話，我們就沒辦法直接計算兩種硬幣正面朝上的概率了，這種情況叫不完全信息。

例如上圖的數據是和完全信息的情況一樣的，區別在于左邊的標簽是問號，不知道是什么硬幣。

這個時候就用到了EM算法。

?

0x23 完全信息下的求解

每次拋硬幣都是獨立的，從二項分布的期望公式$E[\mathbf{X}]=np$可以推導出$p=\frac{E[\mathbf{X}]}{n}=\frac{n_{head}}{n}$。

• 對于A來說，一共拋了三組共三十次，共24次向上6次向下，那么A硬幣朝上的概率是${\hat{\theta }}_A=\frac{24}{30}=0.8$。

• 對于B來說，一共拋了兩組共二十次，共9次向上11次向下，那么B硬幣朝上的概率是${\hat{\theta }}_B=\frac{9}{30}=0.45$。

?

0x24 不完全信息下的初級EM求解

不完全信息情況呢？我們根本不知道每一組的結果是屬于哪種硬幣的，沒辦法用0x24的方法算。這個時候硬幣是否屬于A的隱變量$z_n$是未知的。

（硬幣的情況來說正常用二分法，$ z_n= \begin{cases} 1, if \ coin \ A \ 0, if \ coin \ B \end{cases}$，不過下面使用$z_n$代表數據重新分割的時候屬于A的比例。）

?
?
?

那怎么辦？

想一下就發現，一組拋多次，不同硬幣的拋出不同結果的概率是相當不同的。比如說：

• 一個$\theta =0.3$的硬幣拋出4H6T的概率是$P(4H6T|\theta =0.3)=(\genfrac{}{}{0ex}{}{10}{4})0.3^4(1?0.3{)}^6=\left(\genfrac{}{}{0ex}{}{10}{4}\right)0.0009529569$
• 而$\theta =0.4$的硬幣拋出4H6T的概率是$P(4H6T|\theta =0.4)=(\genfrac{}{}{0ex}{}{10}{4})0.4^4(1?0.4{)}^6=\left(\genfrac{}{}{0ex}{}{10}{4}\right)0.0011943936$

也就是說，倒過來說，看到4H6T的結果的時候，這個硬幣本身朝上的概率更有可能是$\theta =0.4$而不是$\theta =0.3$。

（注意有些文章里面的概率函數式子的寫法用到了分號，$P(\theta ;X)$，意思這是個以X為輸入以$\theta$為變量的函數。為了方便，本文不使用；符號。）

所以說，已知${\theta }_A$和${\theta }_B$的話，我們可以通過觀察拋出來的結果來推測原來硬幣究竟是屬于A還是B的！（這個做法叫做最大似然估計）

可是我們現在不知道${\theta }_A$和${\theta }_B$怎么辦呢？這不是要求解的參數嗎？

面對這個蛋生雞還是雞生蛋的cul-de-sac（死胡同），我們的做法是：先蒙一個！然后再不停互相更新修改。

?
?
?

具體步驟

（1）先給${\theta }_A$和${\theta }_B$隨便賦值。
比如${\theta }_A^{(0)}=0.60$, ${\theta }_B^{(0)}=0.50$。

（2）然后算出

• A硬幣拋出第一組的似然函數是$P(5H5T|{\theta }_A^{(0)}=0.6)=(\genfrac{}{}{0ex}{}{10}{5})0.6^5(1?0.6{)}^5=\left(\genfrac{}{}{0ex}{}{10}{5}\right)0.0007962624$。
• B硬幣拋出第一組的似然函數是$P(5H5T|{\theta }_B^{(0)}=0.5)=(\genfrac{}{}{0ex}{}{10}{5})0.5^5(1?0.5{)}^5=\left(\genfrac{}{}{0ex}{}{10}{5}\right)0.0009765625$。

由此可以看到這組比較有可能是屬于A。這個例子先按照比例來把第一組數據劃分給A和B。

• 劃分給A的比例是$z_1=\frac{P(5H5T|{\theta }_A^{(0)})}{P(5H5T|{\theta }_A^{(0)})+P(5H5T|{\theta }_B^{(0)})}\approx 0.45$。
• 同理劃分給B的比例是$1?z_1=\frac{P(5H5T|{\theta }_B^{(0)})}{P(5H5T|{\theta }_A^{(0)})+P(5H5T|{\theta }_B^{(0)})}\approx 0.55$。

對其他組也進行推算，得到$z_2\approx 0.80$，$z_3\approx 0.73$，$z_4\approx 0.35$，$z_5\approx 0.65$。

（3）接下來得到了新的劃分后的數據，可以更新參數了

• 對于A來說，一共有21.3次向上8.6次向下，那么A硬幣朝上的概率是${\hat{\theta }}_A^{(1)}=\frac{21.3}{21.3+8.6}\approx 0.71$。
• 對于B來說，一共有11.7次向上8.4次向下，那么B硬幣朝上的概率是${\hat{\theta }}_B^{(1)}=\frac{11.7}{11.7+8.4}\approx 0.58$。

（4）重復步驟（2）和（3），直到收斂，可以算得第十次循環之后${\hat{\theta }}_A^{(10)}\approx 0.80$，${\hat{\theta }}_B^{(10)}\approx 0.52$。

可以看到這個結果也跟之前完全信息算出來的比較接近。

?
?
?

0x30 EM算法的公式推導

0x31 定義

• m個互相獨立的樣本組成的數據集$\mathbf{X}=({\mathbf{x}}^{(1)},{\mathbf{x}}^{(2)},...{\mathbf{x}}^{(m)})$(這里每個${\mathbf{x}}^{(k)}$對應硬幣例子里面的一組共拋十次的數據，不知道每組屬于哪種)
• 相對應的隱參數$\mathbf{z}=(z^{(1)},z^{(2)},...z^{(m)})$（每組數據屬于哪種硬幣的標記）
• 樣本本身的模型參數$\theta$（硬幣例子就是$\theta =({\theta }_A,{\theta }_B)$)

對應似然函數為

• 觀察到${\mathbf{x}}^{(k)}$的似然函數為$P({\mathbf{x}}^{(k)}|\theta )$（例如硬幣例子的$P(4H6T|\theta =0.3)$）
• 完全信息情況下的似然函數則是$P({\mathbf{x}}^{(k)},z^{(k)}|\theta )$（例如硬幣例子的$P(4H6T,z=1|{\theta }_1=0.3,{\theta }_2=0.4)$。）
  ?
  ?
  ?

0x32 最大似然估計

那么為了求模型參數$\theta$，將$\theta$看成是參數，求解讓各個樣本的似然函數的乘積$ L(\mathbf{\theta})$最大即可。

• 也就是$ \mathbf{\theta}=\mathop{\arg\max}{\mathbf{\theta}} L(\mathbf{\theta})= \mathop{\arg\max}{\mathbf{\theta}} \prod_{k=1}^{{m}P(\mathbf{x}}{(k)}|\mathbf{\theta})= \mathop{\arg\max}{\mathbf{\theta}} \sum{k=1}^{m} \log P(\mathbf{x}^{(k)}|\mathbf{\theta})$，讓$ L(\mathbf{\theta})$對$\mathbf{\theta}$求導為零容易算出$\mathbf{\theta}$
• 如果有隱函數的話則是$ \mathbf{\theta},z=\mathop{\arg\max}{\mathbf{\theta},z} L(\mathbf{\theta},z)= \mathop{\arg\max}{\mathbf{\theta},z} \sum_{k=1}^{m} \log \sum_{z} P(\mathbf{x}^{(k)},z{(k)}|\mathbf{\theta})$，由于包含了$\log \sum_{z}P(\mathbf{x}^{(k)},z{(k)}|\mathbf{\theta})$這個時候求導的計算量就很繁雜了

解決思路是利用Jensen不等式$E[f(x)]\ge f(E(x))$,

將$\log {\sum }_{z}P({\mathbf{x}}^{(k)},z^{(k)}|\theta )$變成$ \sum_{z} \log P(\mathbf{x}^{(k)},z{(k)}|\mathbf{\theta})$，即是下面的(1)到(2)，同時對P乘以了$\frac{q(z^{(k)})}{q(z{(k)})}=1$。

因此，

$\begin{array}{rl}\sum _{k=1}^m\log \sum _{z}P({\mathbf{x}}^{(k)},z^{(k)}|\theta )& =\sum _{k=1}^m\log \sum _{z}q(z^{(k)})\frac{P({\mathbf{x}}^{(k)},z^{(k)}|\theta )}{q(z^{(k)})}(1)\\ & \ge \sum _{k=1}^m\sum _{z}q(z^{(k)})\log \frac{P({\mathbf{x}}^{(k)},z^{(k)}|\theta )}{q(z^{(k)})}(2)\\ & \to \sum _{k=1}^m\sum _{z}q(z^{(k)}|{\mathbf{x}}^{(k)},{\theta }^{old})\log \frac{P({\mathbf{x}}^{(k)},z^{(k)}|\theta )}{q(z^{(k)}|{\mathbf{x}}^{(k)},{\theta }^{old})}(3)\\ & =\sum _{k=1}^m\sum _{z}q(z^{(k)}|{\mathbf{x}}^{(k)},{\theta }^{old})\log P({\mathbf{x}}^{(k)},z^{(k)}|\theta )\\ & ?\sum _{k=1}^m\sum _{z}q(z^{(k)}|{\mathbf{x}}^{(k)},{\theta }^{old})\log q(z^{(k)}|{\mathbf{x}}^{(k)},{\theta }^{old})(4)\\ & =Q(\theta ,{\theta }^{old})+constant(5)\end{array}$

上面出來傳說中的Q輔助函數，讓Q最大化得出新的$\theta$就是所謂的M步。從(2)到(3)步其實是E步。

所以EM算法就是上面推導公式的(3)(4)(5)之間不斷循環直到收斂。

?
?

0x33 意義解讀

?

E步來看,

（以下參考了人人都懂EM算法，略有修改）

(1)右邊乘以了$ \frac{q(z^{(k)})}{q(z{(k)})}=1$，而引進的未知分布q滿足$ \sum_{z} q(z^{(k)})=1$

(2)里面的$ \sum_{z} q(z^{(k)}) \log \frac{P(\mathbf{x}^{(k)},z{(k)}|\mathbf{\theta})}{q(z^{(k)})}$其實是對$ \log \frac{P(\mathbf{x}^{(k)},z{(k)}|\mathbf{\theta})}{q(z^{(k)})}$求加權平均，也就是求它的數學期望（Expectation）：$ E(\log \frac{P(\mathbf{x}^{(k)},z{(k)}|\mathbf{\theta})}{q(z^{(k)})})$，這也是E步的名字來源。

為了讓(2)能夠取等號，也就是讓$ L(\mathbf{\theta},z)$取一個下限，Jensen不等式告訴我們上面的數學期望里面的變量需要是一個常數，即$ \log \frac{P(\mathbf{x}^{(k)},z{(k)}|\mathbf{\theta})}{q(z^{(k)})}=c$

去掉log之后有$P({\mathbf{x}}^{(k)},z^{(k)}|\theta )=cq(z^{(k)})$

累加后${\sum }_{z}P({\mathbf{x}}^{(k)},z^{(k)}|\theta )=c{\sum }_{z}q(z^{(k)})=c$

所以即$ q(z^{{(k)})=\frac{P(\mathbf{x}}{(k)},z^{{(k)}|\mathbf{\theta})}{c}=\frac{P(\mathbf{x}}{(k)},z^{{(k)}|\mathbf{\theta})}{\sum_{z}P(\mathbf{x}}{(k)},z^{{(k)}|\mathbf{\theta})}=\frac{P(\mathbf{x}}{(k)},z^{{(k)}|\mathbf{\theta})}{P(\mathbf{x}}{(k)}|\mathbf{\theta})}=P(z^{{(k)}|\mathbf{x}}{(k)},\mathbf{\theta})$

q也就是已知$\theta$和${\mathbf{x}}^{(k)}$情況下求隱變量${\mathbf{z}}^{(k)}$的分布，也就是隱變量的后驗概率。然后我們才能繼續算下面M步需要用到的$Q(\theta ,{\theta }^{old})$

$\theta$不是未知數么？說得好，所以(3)里面代入了上次迭代算出的模型參數${\theta }^{old}$。

$Q(\theta ,{\theta }^{old})=E_{\mathbf{z}|\mathbf{X},{\theta }^{old}}(\log L(\mathbf{X},\mathbf{Z}|\theta ))$

以上是E步。

?

M步來看，

最大化$Q(\theta ,{\theta }^{old})$更新$\theta$。（(5)右邊的constant可以忽略，不影響最大化似然函數的操作）

也就是$\theta ={\mathrm{argmax}}_{\theta }Q(\theta ,{\theta }^{old})={\mathrm{argmax}}_{\theta }{\sum }_{k=1}^m{\sum }_{z}q(z^{(k)}|{\mathbf{x}}^{(k)},{\theta }^{old})\log P({\mathbf{x}}^{(k)},z^{(k)}|\theta )$

以上是M步。

?

另一個角度來看，

$\begin{array}{rl}\sum _{k=1}^m\log P({\mathbf{x}}^{(k)}|\theta )& =\sum _{k=1}^m\sum _{z}q(z^{(k)})\log \frac{P({\mathbf{x}}^{(k)}|\theta )P(z^{(k)}|{\mathbf{x}}^{(k)},\theta )}{P(z^{(k)}|{\mathbf{x}}^{(k)},\theta )}(6)\\ & =\sum _{k=1}^m\sum _{z}q(z^{(k)})\log \frac{P({\mathbf{x}}^{(k)},z^{(k)}|\theta )}{P(z^{(k)}|{\mathbf{x}}^{(k)},\theta )}(7)\\ & =\sum _{k=1}^m\sum _{z}q(z^{(k)})\{\log \frac{P({\mathbf{x}}^{(k)},z^{(k)}|\theta )}{q(z^{(k)})}?\log \frac{P(z^{(k)}|{\mathbf{x}}^{(k)},\theta )}{q(z^{(k)})}\}(8)\\ & =\sum _{k=1}^m\sum _{z}q(z^{(k)})\log \frac{P({\mathbf{x}}^{(k)},z^{(k)}|\theta )}{q(z^{(k)})}?\sum _{k=1}^m\sum _{z}q(z^{(k)})\log \frac{P(z^{(k)}|{\mathbf{x}}^{(k)},\theta )}{q(z^{(k)})}(9)\\ & =L(\theta ,z)+KL[q(\mathbf{z})||P(\mathbf{z}|\mathbf{x},\theta )](10)\end{array}$

(6)引入了隱函數后，(7)通過條件概率公式變換概率函數，然后就可以得到(10)。這里可以看出，我們是在構造一個隱函數的分布q。因為我們想讓似然函數最大，那就是說(10)第二項的KL散度盡可能小，也就是要讓構造出來的q盡可能和真實的隱函數分布接近，這時候$q(\mathbf{z})=P(\mathbf{z}|\mathbf{x},\theta )$，KL散度為零。

同時可以看出來剛才的(2)步的讓Jensen不等式取等號的操作也是在讓KL散度為零構造下限，也就是讓$q(\mathbf{z})=P(\mathbf{z}|\mathbf{x},\theta )$取q分布的期望（Expectation）當做隱函數分布的估算。

然后利用q分布再對對數似然函數最大化（Maximization）更新$\theta$。

這個也是所謂的九層境界里面的第二層。（EM算法的九層境界：Hinton和Jordan理解的EM算法）

?

0x40 GMM混合高斯分布的例子

有空再更新

其他

參考及延伸

[1] PRML

[2] 知乎: 怎么通俗易懂地解釋EM算法并且舉個例子?:彭一洋的回答有概括性的數學公式

[3] Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm? Nature Biotechnology, 26(8), 897–899. https://doi.org/10.1038/nbt1406

[4] https://ibug.doc.ic.ac.uk/media/uploads/documents/expectation_maximization-1.pdf

[5] 如何感性地理解EM算法？（拋硬幣的詳解）

[6] EMアルゴリズム徹底解説

[7] EM算法的九層境界：?Hinton和Jordan理解的EM算法

[8] 機器學習系列-強填EM算法在理論與工程之間的鴻溝（上）

[9] 機器學習系列-強填EM算法在理論與工程之間的鴻溝（下）

? H?P ?

總結

以上是生活随笔為你收集整理的EM算法从直观到数学理解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。