EM算法

《統計學習方法》第9章 EM算法及其推廣

EM算法是一種迭代算法，用於含有隱含變量（hidden variable）的機率模型參數的極大似然估計，或極大後驗機率估計。

EM算法的每次迭代由兩步組成：

• E步，求指望（expectation）；
• M步，求極大（maximization）。

因此這一算法成爲指望極大算法（expectation maximization algorithm），簡稱EM算法。

1. EM算法的引入

機率模型有時既含有觀測變量（observable variable），又含有隱變量或潛在變量（latent variable）。若是機率模型的變量都是觀測變量，那麼給定數據，能夠直接用極大似然估計法，或貝葉斯估計發估計模型參數。可是，當模型含有隱變量時，就不能簡單地使用這些估計方法。EM算法就是含有隱變量的機率模型參數的極大似然估計法，或極大後驗機率估計法。咱們在這裏只討論極大似然估計，極大後驗機率估計與其相似。

1.1 EM算法

首先介紹一個使用EM算法的例子。

三硬幣模型 假設有3枚硬幣，分別記做A，B，C。這些硬幣正面出現的機率分別是$ \pi $、$ p $、$ q $。進行以下拋硬幣試驗：先擲硬幣A，根據其結果選出硬幣B或硬幣C，正面選硬幣B，反面選硬幣C；而後擲選出的硬幣，擲硬幣的結果，出現正面記做1，出現反面記做0；獨立地重複n次試驗（這裏，n=10），觀測結果以下：

1, 1, 0, 1, 0, 0, 1, 0, 1, 1

假設只能觀測到擲硬幣的結果，不能觀測擲硬幣的過程。問如何估計三硬幣正面出現的機率，即三硬幣模型的參數。

解 三硬幣模型能夠寫做

$ P(y|\theta)=\sum_{z}{P(y,z|\theta)}=\sum_{z}{P(z|\theta)P(y|z,\theta)}=\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y} $

這裏，隨機變量y是觀測變量，表示一次試驗觀測的結果是1或0；隨機變量z是隱變量，表示未觀測到的擲硬幣A的結果；$ \theta=(\pi,p,q) $是模型參數。這一模型是以上數據的生成模型。

將觀測數據表示爲$ Y=(Y_1,Y_2,…,Y_n)^T $，未觀測數據表示爲$ Z=(Z_1,Z_2,…,Z_n)^T $，則觀測數據的似然函數爲：

$ P(Y|\theta)=\sum_{Z}{P(Z|\theta)P(Y|Z,\theta)} $

即

$ P(Y|\theta)=\prod_{j=1}^{n}{[\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]} $

考慮求模型參數$ \theta=(\pi,p,q) $的極大似然估計，即

$ \hat{\theta}=\arg\max_{\theta}{\log P(Y|\theta)} $

這個問題只能經過迭代的方法求解。

EM算法首先選取參數的初值，記做$ \theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)}) $，而後經過下面的步驟迭代計算參數的估計值，直至收斂爲止。

第i次迭代參數的估計值爲$ \theta^{(i)}=(\pi^{(i)},p^{(0)},q^{(i)}) $，EM算法的第i+1次迭代以下：

 

E步：計算在模型參數$ \pi^{(i)} $，$ p^{(i)} $，$ q^{(i)} $下觀測數據$ y_j $來自擲硬幣B的機率

$ \mu^{(i+1)}= \frac{\pi^{i}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}{\pi^{i}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}+(1-\pi^{i})(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}} $

M步：計算模型參數的新估計值

$ \pi^{(i+1)}=\frac{1}{n}\sum_{j=1}^{n}{\mu_j^{(i+1)}} $

$ p^{(i+1)}=\frac{\sum_{j=1}^{n}{\mu_j^{(i+1)}y_j}}{\sum_{j=1}^{n}{\mu_j^{(i+1)}}} $

$ q^{(i+1)}=\frac{\sum_{j=1}^{n}{(1-\mu_j^{(i+1)})y_j}}{\sum_{j=1}^{n}{(1-\mu_j^{(i+1)})}} $

 

假設模型參數的初值取爲

$ \pi^{(0)}=0.5 $, $ p^{(0)}=0.5 $, $ q^{(0)}=0.5 $

利用迭代公式，獲得模型參數$\theta$的極大似然估計：

$ \hat{\pi}=0.5 $ , $ \hat{p}=0.6 $, $ \hat{q}=0.6 $

 

通常地，用Y表示觀測隨機變量的數據，Z表示隱隨機變量的數據。Y和Z連在一塊兒稱爲徹底數據（complete-data），觀測數據Y又稱爲不徹底數據（incomplete-data）。假設給定觀測數據Y，其機率分佈是$ P(Y|\theta) $，其中$ \theta $是須要估計的模型參數，那麼不徹底數據Y的似然函數是$ P(Y|\theta) $，對數似然函數$ L(\theta)=\log P(Y|\theta) $；假設Y和Z的聯合機率分佈是$ P(Y,Z|\theta) $，那麼徹底數據的對數似然函數是$ \log P(Y,Z|\theta) $。

EM算法經過迭代求$ L(\theta)=\log P(Y|\theta) $的極大似然估計。每次迭代包括兩步：E步，求指望；M步，求極大化。下面來介紹EM算法。