EM算法

一般地，用$Y$表示观测随机变量的数据，$Z$表示隐随机变的数据。$Y$和$Z$连在一起称为完全数据。假设给定观测数据$Y$，其概率分布是$P(Y|\theta )$，其中$\theta$是需要估计地模型参数。

EM算法通过迭代求$L(\theta )=logP(Y|\theta )$的极大似然估计。每次迭代包含两步：$E$步，求期望；$M$步，求极大化。

1. 算法流程

$输入：观测变量数据Y，隐变量数据Z，联合分布P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$

$输出：模型参数\theta$

1. 选择参数的初值${\theta }^{(0)}$，开始迭代
2. $E$步：记${\theta }^{(i)}$为第$i$次迭代参数的估计值，在第$i+1$次迭代的$E$步，计算
   $$Q(\theta ,{\theta }^{(i)})=E_z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
   其中$P(Z|Y,{\theta }^{(i)})$是在给定观测数据$Y$和当前参数估计${\theta }^{(i)}$下隐变量数据$Z$的条件概率分布
3. $M$步：求使$Q(\theta ,{\theta }^{(i)})$的极大化的$\theta$，确定第$i+1$次迭代的参数估计值${\theta }^{(i+1)}$
   $${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
4. 重复第2步和第3步，直到收敛

函数$Q(\theta ,{\theta }^{(i)})$是EM算法的核心，称为Q函数

完全数据的对数似然函数$logP(Y,Z|\theta )$关于在给定观测数据$Y$和当前参数${\theta }^{(i)}$下对未观测数据$Z$的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望称为Q函数

2. EM算法的导出

对于含有隐变量的概率模型，目标是极大化观测数据$Y$关于参数$\theta$的对数似然函数，即极大化
$$L(\theta )=logP(Y|\theta )=log\sum _{Z}P(Y,Z|\theta )=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )$$
$EM$算法是通过迭代逐步近似极大化$L(\theta )$。假设在第$i$次迭代后$\theta$的估计值是${\theta }^{(i)}$，我们希望新估计值$\theta$能使$L(\theta )$增加，即$L(\theta )>L({\theta }^{(i)})$，并逐步到达极大值。为此，考虑两者的差：

$$L(\theta )-L({\theta }^{(i)})=log\sum _{Z}P(Y|Z,\theta )P(Z|\theta )-logP(Y|{\theta }^{(i)})$$

利用Jensen不等式：

$$log\sum _j{\lambda }_j{y}_i\ge \sum _j{\lambda }_{j}log{y}_i,{\lambda }_j\ge 0,\sum _j{\lambda }_j=1$$

得到下界：
$$\begin{gathered} L(\theta )-L({\theta }^{(i)})=log(\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})})-logP(Y|{\theta }^{(i)}) \\ \ge \sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|{\theta }^{(i)}) \\ =\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(Y|{\theta }^{(i)})} \end{gathered}$$

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

则

$$L(\theta )\ge B(\theta ,{\theta }^{(i)})$$

因此，任何可以使$B(\theta ,{\theta }^{(i)})$增大的$\theta$，也可以使$L(\theta )$增大。为了使$L(\theta )$有尽可能大的增大，选择${\theta }^{(i+1)}$使$B(\theta ,{\theta }^{(i)})$达到极大，即
$${\theta }^{(i+1)}=arg\underset{\theta }{\max }B(\theta ,{\theta }^{(i)})$$
$$\begin{gathered} {\theta }^{(i+1)}=arg\underset{\theta }{\max }(L({\theta }^{(i)})+\sum _{Z}P(Z|Y,{\theta }^{(i)})log\frac{P(Y|Z,\theta )P(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)logP(Y,Z|\theta ))=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)}) \end{gathered}$$

3. EM算法在高斯混合模型学习中的应用

3.1 高斯混合模型

高斯混合模型是指具有如下形式的概率分布模型：
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$
其中，${\alpha }_k$是系数，${\alpha }_k\ge 0,{\sum }_{k=1}^K{\alpha }_k=1;\varphi (y|{\theta }_k)$是高斯密度函数，${\theta }_k=({\mu }_k,{\sigma }_k)$
$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}exp(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2})$$
称为第$k$个模型

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

3.2.1 问题描述

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

3.2.2 明确隐变量，写出完全数据的对数似然函数

观测数据$y_j,j=1,2,...,N$是先以概率${\alpha }_k$选择第$k$个高斯分布模型$\varphi (y|{\theta }_k)$，然后通过这个高斯分布模型生成观测数据$y_j$

观测数据$y_j$是已知的，反映观测数据$y_j$来自第$k$个高斯分布模型是未知的，以隐变量${\gamma }_{jk}$表示，定义如下：
$$\begin{gathered} {\gamma }_{jk}=\{\begin{array}{ll}1,& \text{第j个观测数据来自第k个高斯分布模型}\\ 0,& \text{other}\end{array} \\ j=1,2,...,N;k=1,2,...,K \end{gathered}$$
有了观测数据$y_j$和未观测数据${\gamma }_{jk}$，那么完全数据是
$$(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK})$$
完全数据的似然函数：
$$\begin{gathered} P(y,\gamma |\theta )=\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK}|\theta ) \\ =\prod _{j=1}^N\prod _{k=1}^K[{\alpha }_k\varphi (y_j|\theta ){]}^{{\gamma }_{jk}} \\ =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\varphi (y_j|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-{\mu }_k{)}^2}{2{\sigma }_k^2}){]}^{{\gamma }_{jk}} \end{gathered}$$
其中，$n_k={\sum }_{j=1}^N{\gamma }_{jk},{\sum }_{k=1}^K{n}_k=N$

那么完全数据的对数似然函数为：
logP(y,γ∣θ)=∑k=1K{nklogαk+∑j=1Nγjk[log(12π)−logσk−12σk2(yj−μk)2]}logP(y,\gamma|\theta)=\sum_{k=1}^K\{n_klog\alpha_k+\sum_{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]\}logP(y,γ∣θ)=k=1∑K​{nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]}

3.2.3 EM算法的E步：确定Q函数

Q(θ,θ(i))=Eγ[logP(y,γ∣θ)∣y,θ(i)]=Eγ{∑k=1K{nklogαk+∑j=1Nγjk[log(12π)−logσk−12σk2(yj−μk)2]}}=∑k=1K{nklogαk+∑j=1N(Eγjk)[log(12π)−logσk−12σk2(yj−μk)2]}Q(\theta,\theta^{(i)})=E_\gamma[logP(y,\gamma|\theta)|y,\theta^{(i)}]\\=E_\gamma\{\sum_{k=1}^K\{n_klog\alpha_k+\sum_{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\{n_klog\alpha_k+\sum_{j=1}^N(E\gamma_{jk})[log(\frac {1} {\sqrt {2\pi}})-log\sigma_k-\frac{1}{2\sigma_k^2}(y_j-\mu_k)^2]\}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​{nk​logαk​+j=1∑N​(Eγjk​)[log(2π​1​)−logσk​−2σk2​1​(yj​−μk​)2]}

其中
$$\begin{gathered} E({\gamma }_{rk}|y,\theta )=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_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)},j=1,...,N;k=1,..,K \end{gathered}$$

$E({\gamma }_{jk})$表示在当前模型参数下第$j$个观测数据来自第$k$个高斯分布模型的概率，称为第$k$个高斯分布模型对观测数据$y_j$的响应度。

3.2.4 确定EM算法的M步

$${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
分别对${\mu }_k,{\sigma }_k^2$求偏导并令其为0
$$\begin{gathered} \frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\mu }_k}=\sum _{j=1}^{N}E({\gamma }_{jk})\frac{y_i-{\mu }_k}{{\sigma }_k^2}=0 \\ \sum _{j=1}^{N}E({\gamma }_{jk})(y_j-{\mu }_k)=0 \\ {\hat{\mu }}_k=\frac{{\sum }_{j=1}^{N}E({\gamma }_{jk})y_j}{{\sum }_{j=1}^{N}E({\gamma }_{jk})} \end{gathered}$$

---

$$\begin{gathered} \frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\sigma }_k^2}=\sum _{j=1}^{N}E({\gamma }_{jk})[-\frac{1}{2{\sigma }_k^2}+\frac{(y_j-{\mu }_k{)}^2}{2({\sigma }_k^2{)}^2}]=0 \\ \sum _{j=1}^{N}E({\gamma }_{jk})[(y_j-{\mu }_k{)}^2-{\sigma }_k^2]=0 \\ {\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^{N}E({\gamma }_{jk})(y_j-{\mu }_k{)}^2}{{\sum }_{j=1}^{N}E({\gamma }_{jk})} \end{gathered}$$

---

在${\sum }_{k=1}^K{\alpha }_k=1$的条件下对${\alpha }_k$求偏导并令其为0。采用拉格朗日乘子法，有：
$$L(\theta )=Q(\theta ,{\theta }^{(i)})+\lambda (\sum _{k=1}^K{\alpha }_k-1)$$

$$\begin{gathered} \frac{\partial L(\theta )}{\partial {\alpha }_k}=\frac{n_k}{{\alpha }_k}+\lambda =0⟺{\hat{\alpha }}_k=-\frac{n_k}{\lambda } \\ \sum _{k=1}^K(n_k+{\alpha }_k\lambda )=0 \\ N+\lambda =0⟺\lambda =-N \\ {\hat{\alpha }}_k=\frac{n_k}{N}=\frac{{\sum }_{j=1}^{N}E({\gamma }_{jk})}{N} \end{gathered}$$

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

$输入：观测数据y_1,y_2,...,y_N，高斯混合模型$

$输出：高斯混合模型参数$

1. 取参数的初始值开始迭代
2. $E$步：依据当前模型参数，计算高斯分布模型$k$对观测数据$y_j$的响应度
   $$E({\gamma }_{rk}|y,\theta )=\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)}$$
3. $M$步：计算新一轮迭代的模型参数
   $$\begin{gathered} {\hat{\mu }}_k=\frac{{\sum }_{j=1}^{N}E({\gamma }_{jk})y_j}{{\sum }_{j=1}^{N}E({\gamma }_{jk})} \\ {\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^{N}E({\gamma }_{jk})(y_j-{\mu }_k{)}^2}{{\sum }_{j=1}^{N}E({\gamma }_{jk})} \\ {\hat{\alpha }}_k=\frac{{\sum }_{j=1}^{N}E({\gamma }_{jk})}{N} \end{gathered}$$
4. 重复第2步和第3步，直到收敛