机器学习笔记五——EM算法与高斯混合模型

一、EM算法

概率模型有时既含有观测变量，又含有隐变量(latent variable)。如果只含观测变量，那么直接用极大似然估计法估计模型参数即可；但当模型含有隐变量时，就需要采用EM算法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。

1.1 EM算法

用$Y$表示观测随机变量的数据，$Z$表示隐随机变量的数据，$Y$和$Z$连在一起称为完全数据。假设给定观测数据$Y$，其概率分布是$P(Y|\theta )$，其中$\theta$是需要估计的模型参数，$Y$和$Z$的联合概率分布是$P(Y,Z|\theta )$，那么完全数据的对数似然函数是$\log P(Y,Z|\theta )$，$P(Z|Y,\theta )$是给定观测数据$Y$和参数$\theta$下隐变量$Z$的条件概率分布，EM算法通过迭代求解$L(\theta )=\log P(Y|\theta )$的极大似然估计：

$\bullet$选择参数的初始值${\theta }^{(0)}$
$\bullet E$步：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次，计算$$Q(\theta ,{\theta }^{(i)})=E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\text{\hspace{1em}(1)}$$
其中，$P(Z|Y,{\theta }^{(i)})$是在给定观测数据$Y$和当前的参数估计${\theta }^{(i)}$下隐变量数据$Z$的条件概率分布，$Q(\theta ,{\theta }^{(i)})$的第1个变元表示要极大化的参数，第二个变元表示参数的当前估计值，每次迭代实际在求$Q(\theta ,{\theta }^{(i)})$及其极大。

$\bullet M$步：求使$Q(\theta ,{\theta }^{(i)})$极大化的$\theta$，确定第$i+1$步迭代的参数的估计值${\theta }^{(i+1)}$$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(2)}$$

$\bullet$重复$E$步和$M$步，直至收敛

$Q$函数：
完全数据的对数似然函数$\log P(Y,Z|\theta )$ 关于在给定观测数据$Y$和当前参数${\theta }^{(i)}$下 对未观测数据$Z$的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望 称为$Q$函数，即$$Q(\theta ,{\theta }^{(i)})=E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]\text{\hspace{1em}(3)}$$

1.2 EM算法的导出

为什么EM算法能近似实现对观测数据的极大似然估计呢？

对于含有隐变量的概率模型，目标是极大化观测数据$Y$关于参数$\theta$的对数似然函数，即极大化$$L(\theta )=\log P(Y|\theta )=\log \sum _{Z}P(Y,Z|\theta )=\log (\sum _{Z}P(Z|\theta )P(Y|Z,\theta ))\text{\hspace{1em}(4)}$$
注意到$P(Y,Z|\theta )=P(Z|\theta )P(Y|Z,\theta )$

极大化$(4)$式的主要困难是含有未观测数据并有包含和（或积分）的对数。

EM算法是通过逐步迭代近似极大化$L(\theta )$的。假设${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，我们希望新的估计值$\theta$能使$L(\theta )$增加，即$L(\theta )>L({\theta }^{(i)})$，为此，考虑它们的差$$L(\theta )-L({\theta }^{(i)})=\log (\sum _{Z}P(Z|\theta )P(Y|Z,\theta ))-\log P(Y|{\theta }^{(i)})\text{\hspace{1em}(5)}$$
利用Jensen不等式[1]得到其下界$$L(\theta )-L({\theta }^{(i)})=\log (\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})})-\log P(Y|{\theta }^{(i)})$$$$\ge \sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})}-\log P(Y|{\theta }^{(i)})$$$$=\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}\text{\hspace{1em}(6)}$$

令$B(\theta ,{\theta }^{(i)})=L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})}$，则$L(\theta )\ge B(\theta ,{\theta }^{(i)})$

即函数$B(\theta ,{\theta }^{(i)})$是$L(\theta )$的一个下界，为了使$L(\theta )$尽可能增大，可选择${\theta }^{(i+1)}$使$B(\theta ,{\theta }^{(i)})$尽可能增大，即$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }B(\theta ,{\theta }^{(i)})\text{\hspace{1em}(7)}$$

所以$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }(L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})})$$$$=\arg \underset{\theta }{\max }(\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \frac{P(Z|\theta )P(Y|Z,\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})})$$$$=\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})\log (P(Y|Z,\theta )P(Z|\theta ))$$$$=\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^{(i)})\log (P(Y,Z|\theta ))$$$$=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(8)}$$

EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

1.3 EM算法的收敛

设${\theta }^{(i)}(i=1,2,\cdots )$为EM算法得到参数估计序列，$P(Y|{\theta }^{(i)})$为对应的似然函数序列，则$P(Y|{\theta }^{(i)})$是单调递增的，即$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})\text{\hspace{1em}(9)}$$

EM算法不能保证找到全局最优值，因此初始值的选择会对最后的结果产生影响，常用的办法是选取几个不同的初始值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

二、高斯混合模型

2.1 高斯混合模型

高斯混合模型是指具有如下形式的概率分布模型$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)\text{\hspace{1em}(10)}$$
其中，${\alpha }_k$是系数，${\alpha }_k\ge 0,\sum _{k=1}^K{\alpha }_k=1$，$\varphi (y|{\theta }_k)$是高斯分布密度函数，${\theta }_k=({\mu }_k,{\sigma }_k^2)$ ϕ(y∣θk)=12πσkexp⁡{−(y−μk)22σk2}(11)\phi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}\exp\{-\frac{(y-\mu_k)^2}{2\sigma_k^2}\}\tag{11}ϕ(y∣θk​)=2π​σk​1​exp{−2σk2​(y−μk​)2​}(11)

2.2 高斯混合模型参数估计的EM算法

假设观测数据$y_1,\cdots ,y_N$由高斯混合模型$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$生成，其中，$\theta =({\alpha }_1,\cdots ,{\alpha }_K;{\theta }_1,\cdots ,{\theta }_K)$，我们用EM算法来估计参数$\theta$

定义隐变量$$\begin{gathered} {\gamma }_{jk}=\{\begin{array}{ll}1& 第j个观测数据来自第k个分模型\\ 0& 否则\end{array} \\ \text{\hspace{1em}(12)} \end{gathered}$$
其中$j=1,2\cdots ,N;k=1,2\cdots ,K$

从而有完全数据$(y_j,{\gamma }_{j1},\cdots ,{\gamma }_{jK})$，$j=1,2\cdots ,K$，所以可以写出完全数据的似然函数P(y,γ∣θ)=∏j=1NP(yj,γj1,⋯ ,γjK∣θ)=∏k=1K∏j=1N[αkϕ(yj∣θk)]γjk=∏k=1Kαknk∏j=1N[12πσkexp⁡{−(yj−μk)22σk2}]γjk(12)P(y,\gamma|\theta)=\prod_{j=1}^NP(y_j, \gamma_{j1},\cdots,\gamma_{jK}|\theta)\\ =\prod_{k=1}^K\prod_{j=1}^N[\alpha_k\phi(y_j|\theta_k)]^{\gamma_{jk}}\\ =\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\frac{1}{\sqrt{2\pi}\sigma_k}\exp\{-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}\}]^{\gamma_{jk}}\tag{12}P(y,γ∣θ)=j=1∏N​P(yj​,γj1​,⋯,γjK​∣θ)=k=1∏K​j=1∏N​[αk​ϕ(yj​∣θk​)]γjk​=k=1∏K​αknk​​j=1∏N​[2π​σk​1​exp{−2σk2​(yj​−μk​)2​}]γjk​(12)
其中$n_k=\sum _{j=1}^N{\gamma }_{jk},\sum _{k=1}^K{n}_k=N$
所以对数似然函数为log⁡P(y,γ∣θ)=∑k=1K{nklog⁡αk+∑j=1Nγjk[log⁡(12π)−log⁡σk−12σk2(yj−μk)2]}(13)\log P(y,\gamma|\theta)=\sum_{k=1}^K\{n_k\log\alpha_k+\sum\limits_{j=1}^N\gamma_{jk}[\log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2]\}\tag{13}logP(y,γ∣θ)=k=1∑K​{nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]}(13)

 
$\bullet E-step$
所以Q(θ,θ(i))=E[log⁡P(y,γ∣θ)∣y,θ(i)]=E{∑k=1K{nklog⁡αk+∑j=1Nγjk[log⁡(12π)−log⁡σk−12σk2(yj−μk)2]}=∑k=1K{∑j=1NE(γjk)log⁡αk+∑j=1NE(γjk)[log⁡(12π)−log⁡σk−12σk2(yj−μk)2]}(14)Q(\theta,\theta^{(i)})=E[\log P(y,\gamma|\theta)|y,\theta^{(i)}]=E\{\sum_{k=1}^K\{n_k\log\alpha_k+\sum\limits_{j=1}^N\gamma_{jk}[\log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2]\}\\ =\sum_{k=1}^K\{\sum_{j=1}^NE(\gamma_{jk})\log\alpha_k+\sum_{j=1}^NE(\gamma_{jk})[\log(\frac{1}{\sqrt{2\pi}})-\log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2]\}\tag{14}Q(θ,θ(i))=E[logP(y,γ∣θ)∣y,θ(i)]=E{k=1∑K​{nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]}=k=1∑K​{j=1∑N​E(γjk​)logαk​+j=1∑N​E(γjk​)[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]}(14)
注意到$n_k=\sum _{j=1}^N{\gamma }_{jk}$，所以$E(n_k)=\sum _{j=1}^{N}E({\gamma }_{jk})$

而$$\begin{gathered} E({\gamma }_{jk})=P({\gamma }_{jk}=1|y,\theta ) \\ =\frac{P({\gamma }_{jk}=1,y_j|\theta )}{\sum _{k=1}^{K}P({\gamma }_{jk}=1,y_j|\theta )} \\ =\frac{P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}{\sum _{k=1}^{K}P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )} \\ =\frac{{\alpha }_k\varphi (y_i|{\theta }_k)}{\sum _{k=1}^K{\alpha }_k\varphi (y_i|{\theta }_k)}\text{\hspace{1em}(15)} \end{gathered}$$
其中$j=1,2\cdots ,N;k=1,2\cdots ,K$

所以$E({\gamma }_{jk})$可以看作是在当前模型参数下，第$j$个观测数据来自第$k$个分模型的概率，称为分模型$k$对观测数据$y_j$的响应度。

 
$\bullet M-step$
即${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$。将$Q(\theta ,{\theta }^{(i)})$分别对${\mu }_k,{\sigma }_k^2$分别求偏导并令其为0即可得到${\hat{\mu }}_k,{\hat{\sigma }}_k^2$，在$\sum _{k=1}^K{\alpha }_k=1$的条件下，对${\alpha }_k$求偏导即可得到${\hat{\alpha }}_k$：
$${\hat{\mu }}_k=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{\sum _{j=1}^N{\hat{\gamma }}_{jk}}\text{\hspace{1em}(16)}$$
$${\hat{\sigma }}_k^2=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}(y_j-{\mu }_k{)}^2}{\sum _{j=1}^N{\hat{\gamma }}_{jk}}\text{\hspace{1em}(17)}$$
$${\hat{\alpha }}_k=\frac{n_k}{N}=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{N}\text{\hspace{1em}(18)}$$
其中，$k=1,\cdots ,K$