EM算法

極大似然估計

極大似然估計：(maximum likelihood estimate, MLE) 是一種常用的模型參數估計方法。它假設觀測樣本出現的概率最大，也即樣本聯合概率（也稱似然函數）取得最大值。

為求解方便，對樣本聯合概率取對數似然函數
$$\log L(\theta )=\log \mathbb{P}(X|\theta )=\sum _{i=1}^N\log \mathbb{P}({\mathbf{x}}_i|\theta )$$
優化目標是最大化對數似然函數
$$\hat{\theta }=\arg \underset{\theta }{\max }\sum _{i=1}^N\log \mathbb{P}({\mathbf{x}}_i|\theta )$$

假設瓜田里有兩種類型的西瓜🍉，瓜農隨機抽取了10個西瓜，來了解西瓜的重量分布 $p(x|\theta )$，記錄結果如下：

變量	樣本
西瓜重量 $x$	5.3 , 5.7, 4.7, 4.3, 3.2, 4.9, 4.1, 3.5, 3.8, 1.7
西瓜品種 $z$	1, 1, 1, 1, 2, 2, 2, 2, 2, 2

其中，西瓜的品種 $z$ 是離散分布 $\mathbb{P}(z=k)={\pi }_k$，一般假設兩種類型的西瓜服從均值和方差不同的高斯分布 $N({\mu }_1,{\sigma }_1^2)$和 $N({\mu }_2,{\sigma }_2^2)$。由全概率公式，西瓜重量的概率密度模型
$$p(x;\theta )={\pi }_1\mathcal{N}(x;{\mu }_1,{\sigma }_1^2)+{\pi }_2\mathcal{N}(x;{\mu }_2,{\sigma }_2^2)$$

我們嘗試用極大似然估計求解參數$\theta =({\pi }_1,{\pi }_2,{\mu }_1,{\sigma }_1^2,{\mu }_2,{\sigma }_2^2)$。

優化目標函數
$$\begin{gathered} \underset{\theta }{\max }\sum _{z_i=1}\log {\pi }_1\mathcal{N}(x_i;{\mu }_1,{\sigma }_1^2)+\sum _{z_i=2}\log {\pi }_2\mathcal{N}(x_i;{\mu }_2,{\sigma }_2^2) \\ \text{s.t.?}{\pi }_1+{\pi }_2=1 \end{gathered}$$
使用拉格朗日乘子法容易求得
$$\begin{gathered} {\pi }_1=0.4,{\pi }_2=0.6 \\ {\mu }_1=5,{\sigma }_1^2=0.54^2 \\ {\mu }_2=3.53,{\sigma }_2^2=0.98^2 \end{gathered}$$

最終得到

$$p(x)=0.4\times \mathcal{N}(x;5,0.54^2)+0.6\times \mathcal{N}(x;3.53,0.98^2)$$

但是，實際中如果瓜農無法辯識標記西瓜的品種，此時概率分布函數變為
$$p(x;\theta )=\pi \mathcal{N}(x;{\mu }_1,{\sigma }_1^2)+(1?\pi )\mathcal{N}(x;{\mu }_2,{\sigma }_2^2)$$

其中品種$z$ 成為隱藏變量。對數似然函數變為
$$\log L(\theta )=\sum _i\log (\pi \mathcal{N}(x_i;{\mu }_1,{\sigma }_1^2)+(1?\pi )\mathcal{N}(x_i;{\mu }_2,{\sigma }_2^2))$$
其中參數 $\theta =(\pi ,{\mu }_1,{\sigma }_1^2,{\mu }_2,{\sigma }_2^2)$。上式中存在"和的對數"，若直接求導將會變得很麻煩。下節我們將會介紹EM算法來解決此類問題。

基本思想

概率模型有時既含有觀測變量 (observable variable)，又含有隱變量 (latent variable)。EM（Expectation-Maximization，期望最大算法）是一種迭代算法，用于含有隱變量的概率模型的極大似然估計或極大后驗估計，是數據挖掘的十大經典算法之一。

假設現有一批獨立同分布的樣本
X = { x 1 , x 2 , ? , x N } X=\{x_1,x_2,\cdots,x_N\} X={x1?,x2?,?,xN?}
它們是由某個含有隱變量的概率分布 $p(x,z|\theta )$ 生成。設樣本對應的隱變量數據
Z = { z 1 , z 2 , ? , z N } Z=\{z_1,z_2,\cdots,z_N\} Z={z1?,z2?,?,zN?}

對于一個含有隱變量 $Z$ 的概率模型，一般將 $(X,Z)$ 稱為完全數據 (complete-data)，而觀測數據 $X$ 為不完全數據(incomplete-data)。

假設觀測數據 $X$ 概率密度函數是$p(X|\theta )$，其中$\theta$是需要估計的模型參數，現嘗試用極大似然估計法估計此概率分布的參數。為了便于討論，此處假設 $z$ 為連續型隨機變量，則對數似然函數為
$$\log L(\theta )=\sum _{i=1}^N\log p(x_i|\theta )=\sum _{i=1}^N\log {\int }_{z_i}p(x_i,z_i|\theta )\mathrm{d}{z}_i$$

Suppose you have a probability model with parameters $\theta$.
$p(x|\theta )$ has two names. It can be called the probability of $x$ (given $\theta$), or the likelihood of $\theta$ (given that $x$ was observed).

我們的目標是極大化觀測數據 $X$ 關于參數 $\theta$ 的對數似然函數
$$\hat{\theta }=\arg \underset{\theta }{\max }\log L(\theta )$$

顯然，此時 $\log L(\theta )$ 里含有未知的隱變量 $z$ 以及求和項的對數，相比于不含隱變量的對數似然函數，該似然函數的極大值點較難求解，而 EM 算法則給出了一種迭代的方法來完成對 $\log L(\theta )$ 的極大化。

注意：確定好含隱變量的模型后，即確定了聯合概率密度函數 $p(x,z|\theta )$ ，其中$\theta$是需要估計的模型參數。為便于討論，在此有必要說明下其他已知的概率函數。

聯合概率密度函數
$$p(x,z|\theta )=f(x,z;\theta )$$
觀測變量 $x$ 的概率密度函數
$$p(x|\theta )={\int }_{z}f(x,z;\theta )\mathrm{d}z$$
隱變量 $z$ 的概率密度函數
$$p(z|\theta )={\int }_{x}f(x,z;\theta )\mathrm{d}x$$
條件概率密度函數
$$p(x|z,\theta )=\frac{p(x,z|\theta )}{p(z|\theta )}=\frac{f(x,z;\theta )}{{\int }_{x}f(x,z;\theta )\mathrm{d}x}$$
和
$$p(z|x,\theta )=\frac{p(x,z|\theta )}{p(x|\theta )}=\frac{f(x,z;\theta )}{{\int }_{z}f(x,z;\theta )\mathrm{d}z}$$
下面給出兩種推導方法：一種借助 Jensen 不等式；一種使用 KL 散度。

首先使用 Jensen 不等式推導：使用含有隱變量的全概率公式
$$\begin{array}{rl}\log p(x_i|\theta )& =\log {\int }_{z_i}p(x_i,z_i|\theta )\mathrm{d}{z}_i\\ & =\log {\int }_{z_i}{q}_i(z_i)\frac{p(x_i,z_i|\theta )}{q_i(z_i)}\mathrm{d}{z}_i\\ & =\log {\mathbb{E}}_z\left(\frac{p(x_i,z_i|\theta )}{q_i(z_i)}\right)\\ & ?{\mathbb{E}}_z\left(\log \frac{p(x_i,z_i|\theta )}{q_i(z_i)}\right)\\ & ={\int }_{z_i}{q}_i(z_i)\log \frac{p(x_i,z_i|\theta )}{q_i(z_i)}\mathrm{d}{z}_i\end{array}$$
其中 $q_i(z_i)$ 是引入的第$i$個樣本隱變量$z_i$ 的任意概率密度函數（未知函數），其實 $q$ 不管是任意函數，上式都成立。從后續推導得知，當 $q_i(z_i)$ 是 $z_i$ 的概率密度時，方便計算。

所以
$$\log L(\theta )=\sum _{i=1}^N\log p(x_i|\theta )?B(q,\theta )=\sum _{i=1}^N{\int }_{z_i}{q}_i(z_i)\log \frac{p(x_i,z_i|\theta )}{q_i(z_i)}\mathrm{d}{z}_i$$

其中函數 $B$ 為對數似然的下界函數。下界比較好求，所以我們要優化這個下界來使得似然函數最大。

假設第 $t$ 次迭代時 $\theta$ 的估計值是 ${\theta }^{(t)}$，我們希望第 $t+1$ 次迭代時的 $\theta$ 能使 $\log L(\theta )$ 增大，即
$$\log L({\theta }^{(t)})?\log L({\theta }^{(t+1)})$$

可以分為兩步實現：

• 首先，固定$\theta ={\theta }^{(t)}$ ，通過調整 $q$ 函數使得 $B(q^{(t)},\theta )$ 在 ${\theta }^{(t)}$ 處和 $\log L({\theta }^{(t)})$ 相等；
  $$B(q^{(t)},{\theta }^{(t)})=\log L({\theta }^{(t)})$$

• 然后，固定$q$，優化 ${\theta }^{(t+1)}$ 取到下界函數 $B(q^{(t)},\theta )$ 的最大值。
  $${\theta }^{(t+1)}=\arg \underset{\theta }{\max }B(q^{(t)},\theta )$$

所以
$$\log L({\theta }^{(t+1)})?B(q^{(t)},{\theta }^{(t+1)})?B(q^{(t)},{\theta }^{(t)})=\log L({\theta }^{(t)})$$

因此，EM算法也可以看作一種坐標提升算法，首先固定一個值，對另外一個值求極值，不斷重復直到收斂。

接下來，我們開始求解 $q^{(t)}$ 。Jensen不等式中等號成立的條件是自變量是常數，即
$$\frac{p(x_i,z_i|\theta )}{q_i(z_i)}=c$$
由于假設 $q_i(z_i)$是 $z_i$ 的概率密度函數，所以
$$p(x_i|\theta )={\int }_{z_i}p(x_i,z_i|\theta )\mathrm{d}{z}_i={\int }_{z_i}c{q}_i(z_i)\mathrm{d}{z}_i=c$$

于是
$$q_i(z_i)=\frac{p(x_i,z_i|\theta )}{c}=\frac{p(x_i,z_i|\theta )}{p(x_i|\theta )}=p(z_i|x_i,\theta )$$
可以看到，函數 $q_i(z_i)$ 代表第 $i$ 個數據是 $z_i$ 的概率密度，是可以直接計算的。

最終，我們只要初始化或使用上一步已經固定的 ${\theta }^{(t)}$，然后計算
$$\begin{array}{rl}{\theta }^{(t+1)}& =\arg \underset{\theta }{\max }\sum _{i=1}^N{\int }_{z_i}p(z_i|x_i,{\theta }^{(t)})\log \frac{p(x_i,z_i|\theta )}{p(z_i|x_i,{\theta }^{(t)})}\mathrm{d}{z}_i\\ & =\arg \underset{\theta }{\max }\sum _{i=1}^N{\int }_{z_i}p(z_i|x_i,{\theta }^{(t)})\log p(x_i,z_i|\theta )\mathrm{d}{z}_i\\ & =\arg \underset{\theta }{\max }\sum _{i=1}^N{\mathbb{E}}_{z_i|x_i,{\theta }^{(t)}}[\log p(x_i,z_i|\theta )]\\ & =\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(t)})\end{array}$$

接下來使用 KL 散度推導：使用含有隱變量的條件概率
$$\begin{array}{rl}\log p(x_i|\theta )& =\log \frac{p(x_i,z_i|\theta )}{p(z_i|x_i,\theta )}\\ & ={\int }_{z_i}{q}_i(z_i)\log \frac{p(x_i,z_i|\theta )}{p(z_i|x_i,\theta )}?\frac{q_i(z_i)}{q_i(z_i)}\mathrm{d}{z}_i\\ & ={\int }_{z_i}{q}_i(z_i)\log \frac{p(x_i,z_i|\theta )}{q_i(z_i)}\mathrm{d}{z}_i+{\int }_{z_i}{q}_i(z_i)\log \frac{q_i(z_i)}{p(z_i|x_i,\theta )}\mathrm{d}{z}_i\\ & =B(q_i,\theta )+KL(q_i\Vert p_i)\end{array}$$
同樣 $q_i(z_i)$ 是引入的關于 $z_i$ 的任意概率密度函數（未知函數），函數 $B(q_i,\theta )$ 表示對數似然的一個下界，散度 $KL(q_i\Vert p_i)$描述了下界與對數似然的差距。

同樣為了保證
$$\log L({\theta }^{(t)})?\log L({\theta }^{(t+1)})$$

分為兩步實現：

• 首先，固定$\theta ={\theta }^{(t)}$ ，通過調整 $q$ 函數使得 $B(q^{(t)},\theta )$ 在 ${\theta }^{(t)}$ 處和 $\log L({\theta }^{(t)})$ 相等，即 $KL(q_i\Vert p_i)=0$，于是
  $$q_i(z_i)=p(z_i|x_i,{\theta }^{(t)})$$

• 然后，固定$q$，優化 ${\theta }^{(t+1)}$ 取到下界函數 $B(q^{(t)},\theta )$ 的最大值。
  $${\theta }^{(t+1)}=\arg \underset{\theta }{\max }B(q^{(t)},\theta )$$

算法流程

輸入：觀測數據 $X$，聯合分布 $p(x,z;\theta )$，條件分布$P(z|x,\theta )$

輸出：模型參數$\theta$

EM算法通過引入隱含變量，使用極大似然估計（MLE）進行迭代求解參數。每次迭代由兩步組成：

• E-step：求期望 (expectation)。以參數的初始值或上一次迭代的模型參數 ${\theta }^{(t)}$ 來計算隱變量后驗概率 $p(z_i|x_i,{\theta }^{(t)})$ ，并計算期望(expectation)
  $$Q(\theta ,{\theta }^{(t)})=\sum _{i=1}^N{\int }_{z_i}p(z_i|x_i,{\theta }^{(t)})\log p(x_i,z_i|\theta )\mathrm{d}{z}_i$$

• M-step: 求極大 (maximization)，極大化E步中的期望值，來確定 $t+1$ 次迭代的參數估計值
  $${\theta }^{(t+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(t)})$$

依次迭代，直至收斂到局部最優解。

高斯混合模型

基礎模型

高斯混合模型 (Gaussian Mixture Model, GMM) 數據可以看作是從$K$個高斯分布中生成出來的，每個高斯分布稱為一個組件 (Component)。

引入隱變量 z ∈ { 1 , 2 , ? , K } z\in\{1,2,\cdots,K\} z∈{1,2,?,K}，表示對應的樣本 $x$ 屬于哪一個高斯分布，這個變量是一個離散的隨機變量：
$$\begin{gathered} \mathbb{P}(z=k)={\pi }_k \\ \text{s.t.?}\sum _{k=1}^K{\pi }_k=1 \end{gathered}$$
可將 ${\pi }_k$ 視為選擇第 $k$ 高斯分布的先驗概率，而對應的第$k$ 個高斯分布的樣本概率
$$p(x|z=k)=\mathcal{N}(x;{\mu }_k,{\Sigma }_k)$$

于是高斯混合模型
$$p_M(x)=\sum _{k=1}^K{\pi }_k\mathcal{N}(x;{\mu }_k,{\Sigma }_k)$$

其中 $0?{\pi }_k?1$為混合系數(mixing coefficients)。

高斯混合模型的參數估計是EM算法的一個重要應用，隱馬爾科夫模型的非監督學習也是EM算法的一個重要應用。

EM算法

高斯混合模型的極大似然估計
$$\hat{\theta }=\arg \underset{\theta }{\max }\sum _{i=1}^N\log \sum _{k=1}^K{\pi }_k\mathcal{N}(x_i;{\mu }_k,{\Sigma }_k)$$
其中參數 ${\theta }_k=({\pi }_k,{\mu }_k,{\Sigma }_k)$，使用EM算法估計GMM的參數$\theta$。

依照當前模型參數，計算隱變量后驗概率：由貝葉斯公式知道
$$\begin{array}{rl}P(z_i=k|x_i)& =\frac{P(z_i=k)p(x_i|z_i=k)}{p(x_i)}\\ & =\frac{{\pi }_k\mathcal{N}(x_i;{\mu }_k,{\Sigma }_k)}{{\sum }_{k=1}^K{\pi }_k\mathcal{N}(x_i;{\mu }_k,{\Sigma }_k)}\\ & ={\gamma }_{ik}\end{array}$$

令 ${\gamma }_{ik}$表示第$i$個樣本屬于第$k$個高斯分布的概率。

E-step：確定Q函數

$$\begin{array}{rl}Q(\theta ,{\theta }^{(t)})& =\sum _{i=1}^N\sum _{k=1}^{K}p(z_i=k|x_i,{\mu }^{(t)},{\Sigma }^{(t)})\log p(x_i,z_i=k|\mu ,\Sigma )\\ & =\sum _{i=1}^N\sum _{k=1}^K{\gamma }_{ik}\log {\pi }_k\mathcal{N}(x;{\mu }_k,{\Sigma }_k)\\ & =\sum _{i=1}^N\sum _{k=1}^K{\gamma }_{ik}(\log {\pi }_k+\log \mathcal{N}(x;{\mu }_k,{\Sigma }_k))\end{array}$$

M-step：求Q函數的極大值

上面已獲得的$Q(\theta ,{\theta }^{(t)})$分別對 ${\mu }_k,{\Sigma }_k$求導并設為0。得到
$$\begin{gathered} {\mu }_k^{(t+1)}=\frac{{\sum }_{i=1}^N{\gamma }_{ik}{x}_i}{{\sum }_{i=1}^N{\gamma }_{ik}} \\ {\Sigma }_k^{(t+1)}=\frac{{\sum }_{i=1}^N{\gamma }_{ik}(x_i?{\mu }_k^{(t+1)})(x_i?{\mu }_k^{(t+1)}{)}^T}{{\sum }_{i=1}^N{\gamma }_{ik}} \end{gathered}$$

可以看到第$k$個高斯分布的${\mu }_k,{\Sigma }_k$ 是所有樣本的加權平均，其中每個樣本的權重為該樣本屬于第$k$個高斯分布的后驗概率 ${\gamma }_{ik}$。

對于混合系數 ${\pi }_k$，因為有限制條件，使用拉格朗日乘子法可求得
$${\pi }_k^{(t+1)}=\frac{1}{N}\sum _{i=1}^N{\gamma }_{ik}$$

即第$k$個高斯分布的混合系數是屬于$k$的樣本的平均后驗概率，由此運用EM算法能大大簡化高斯混合模型的參數估計過程，在中間步只需計算${\gamma }_{ik}$就行了。

高斯混合模型的算法流程如下圖所示：

高斯混合聚類

高斯混合聚類假設每個類簇中的樣本都服從一個多維高斯分布，那么空間中的樣本可以看作由$K$個多維高斯分布混合而成。

引入隱變量$z$ 標記簇類別，這樣就可以使用高斯混合模型
$$p_M(x)=\sum _{k=1}^K{\pi }_k\mathcal{N}(x;{\mu }_k,{\Sigma }_k)$$

使用EM算法迭代求解。

相比于K-means更具一般性，能形成各種不同大小和形狀的簇。K-means可視為高斯混合聚類中每個樣本僅指派給一個混合成分的特例
$${\gamma }_{ik}=\{\begin{array}{ll}1,& \text{if?}k=\arg {\min }_k\Vert x_i?{\mu }_k{\Vert }^2\\ 0,& \text{otherwise}\end{array}$$
且各混合成分協方差相等，均為對角矩陣 ${\sigma }^{2}I$。

附錄

Jensen 不等式

若 $f$ 是凸函數(convex function)，對任意的 $\lambda \in [0,1]$，下式恒成立
$$f(\lambda x_1+(1?\lambda )x_2)?\lambda f(x_1)+(1?\lambda )f(x_2)$$
Jensen’s inequality就是上式的推廣，設 $f(x)$ 為凸函數，${\lambda }_i\in [0,1],{\sum }_i{\lambda }_i=1$ ，則
$$f(\sum _i{\lambda }_i{x}_i)?\sum _i{\lambda }_{i}f(x_i)$$
若將 ${\lambda }_i$ 視為一個概率分布，則可表示為期望值的形式
$$f(\mathbb{E}[x])?\mathbb{E}[f(x)]$$

顯然，如果 $f$ 是凹函數(concave function)，則將不等號反向。