EM算法和高斯混合聚类

EM算法引言

在现实应用中，概率模型有时既含有观测变量(observable variable)，又含有不能被观测到的变量，该变量称为隐变量(latent variable)。如果给定数据全都是观测变量，那么可以使用最大似然估计求解模型参数，但是在含有隐变量的情况下无法求解。EM算法就是用于求解在训练样本具有隐变量的情况下概率模型参数的最大似然估计。EM算法是对两种未知参数（隐变量分布和模型参数）交替优化的迭代算法，相当于在优化过程中，保持某一个变量不变，去优化另一个变量。首先给定模型参数，EM算法的每次迭代分为两步：E步，求隐变量的期望（expectation）；M步，求模型参数的极大似然估计（maximization），所以这一算法称为期望极大算法（expectation maximization algorithm）。
实例：

从上面的例子中可以看出，参数的最大似然估计为$P(y|\theta )$，此时由于观测值$y$与$\theta$无直接关系，需要加入隐变量，即为：
$$\begin{gathered} 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} \end{gathered}$$
$P(y|\theta )={\sum }_{z}P(y,z|\theta )$相当于全概公式。

将观测数据表示为$Y=(Y_1,Y_2,...,Y_n)$,隐藏数据表表示为$Z=(Z_1，Z_2,...,Z_N)$,则观测数据的似然函数为：
$$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}]$$
求模型参数$\hat{\theta }=\arg {\max }_{\theta }\log P(Y|\theta )$
上述问题没有解析解，只有通过迭代的方法求解，EM算法就是求解此问题的迭代算法。

EM算法推导

给定观测数据 Y = { Y ( i ) } i = 1 m Y=\{Y^{(i)} \}_{i=1}^{m} Y={Y(i)}i=1m​,则模型的最大似然估计为：
$$L(\theta )=\sum _{i=1}^m\log P(Y^{(i)}|\theta )=\sum _{i=1}^m\log [\sum _{z}P(Y^{(i)}|z,\theta )P(z|\theta )]$$
$$L(\theta )=\log P(Y|\theta )=\log \sum _{Z}P(Y|Z,\theta )P(Z|\theta )$$
Jensen不等式：

EM算法通过逐步迭代最大化$L(\theta )$，假设第$i$次迭代之后的估计值为${\theta }^{(i)}$,我们希望新估计的值$\theta$能使$L(\theta )$增大逐步达到极大值。考虑两者的差值：
$$\begin{gathered} L(\theta )-L({\theta }^{(i)})=\log \sum _{Z}P(Y|Z,\theta )P(Z|\theta )-\log P(Y|{\theta }^{(i)}) \\ =\log \sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}-\log P(Y|{\theta }^{(i)}) \end{gathered}$$
上式中，${\sum }_{Z}P(Z|Y,{\theta }^{(i)})=1$，左式的log里面的值相当于 $Z|Y,{\theta }^{(i)}$的期望。$$E_{Z|Y,{\theta }^{(i)}}(\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})})=\sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})}$$
由Jensen不等式可得：
$$\log E_{Z|Y,{\theta }^{(i)}}(\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})})\ge E_{Z|Y,{\theta }^{(i)}}\log (\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})})$$
所以可得到：
$$\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)})}]-\log P(Y|{\theta }^{(i)}) \\ =\log \sum _{Z}P(Z|Y,{\theta }^{(i)})\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^{(i)})P(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)})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)})}$$
当$\theta ={\theta }^{(i)}$时，log项为0，此时$L({\theta }^{(i)})=B(\theta ,{\theta }^{(i)})$
可以得到目标函数的下界为：：$$L(\theta )\ge B(\theta ,{\theta }^{(i)})$$
$B(\theta ,{\theta }^{(i)})$越大，可以使得$L(\theta )$的下界越大，因为在该次迭代中我们无法得到$L(\theta )$的最大值，我们只能取下界的最大的值：
$${\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)})\log P(Y,Z|\theta )) \\ =\arg \underset{\theta }{\max }{E}_{Z|Y,{\theta }^{(i)}}(\log P(Y,Z|\theta )) \end{gathered}$$
令：
$$Q(\theta ,{\theta }^{(i)})=E_{Z|Y,{\theta }^{(i)}}(\log P(Y,Z|\theta ))$$
则：
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
该公式相当于EM算法对参数$\theta$的一次迭代求解，通过不断求解下界的极大化值来逼近最大化似然函数。在迭代过程中，$L(\theta )$是不断增大的，但是EM算法并不能保证收敛到全局最优值。

EM算法总流程

1.初始化参数${\theta }^0$
2.E步：求似然函数关于$z$的期望：$Q(\theta ,{\theta }^{(i)})=E_{Z|Y,{\theta }^{(i)}}(\log P(Y,Z|\theta )$
3.M步：更新${\theta }^{(i+1)}=\arg {\max }_{\theta }Q(\theta ,{\theta }^{(i)})$

高斯混合聚类

EM算法作为一种求解思想，可以用于多种模型求解，高斯混合聚类是其典型应用。下图表明，在聚类中，一个高斯分布很难进行聚类，此时需要多个高斯分布进行聚类。

对于$n$维空间中的随机变量$\mathbf{x}$遵从多元高斯分布，其概率密度为：
$$p(\mathbf{x})=\frac{1}{(2\pi {)}^{\frac{n}{2}}|\Sigma {|}^{\frac{1}{2}}}\exp (-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }){\Sigma }^{-1}(\mathbf{x}-\mathbf{\mu }{)}^T)$$
对于存在多个高斯分布的情况下，我们定义高斯混合分布：
$$p(\mathbf{x})=\sum _{k=1}^K{\alpha }_{k}p(\mathbf{x}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)$$
该分布有$K$个高斯分布组成，也就是一共有$K$个簇。${\alpha }_k=P(k)>0$为混合系数，表明选择第$i$个成分的概率，${\sum }_{i=1}^K{\alpha }_k=1$。某样本$\mathbf{x}$对应第$k$个成分的后验概率为：
$$P(k|\mathbf{x})=\frac{p(\mathbf{x}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)p(k)}{p(\mathbf{x})}=\frac{p(\mathbf{x}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)p(k)}{{\sum }_{k=1}^K{\alpha }_{k}p(\mathbf{x}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)}$$
给定数据集 D = { x 1 , x 2 , . . . , x m } D=\{\bm{x}_{1},\bm{x}_{2},...,\bm{x}_{m}\} D={x1​,x2​,...,xm​},此时我们可以通过最大似然估计求解模型参数 { α k , μ k , Σ k } k = 1 K \{\alpha_{k},\bm{\mu}_{k},\bm{\Sigma}_{k}\}_{k=1}^{K} {αk​,μk​,Σk​}k=1K​.
$$p(D)=\prod _{i=1}^{m}p({\mathbf{x}}_i)=\prod _{i=1}^m\sum _{k=1}^K{\alpha }_{k}p({\mathbf{x}}_i|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)$$
取对数似然函数：
$$L(D)=\sum _{i=1}^m\log \sum _{k=1}^K{\alpha }_{k}p({\mathbf{x}}_i|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)$$
对于上述带有$\log \sum$的形式，从经验来说求导很困难，此时可以EM算法求解。
首先明确隐变量
观测数据是已知的，但是数据来源于哪个成分是未知的。对于观测数据${\mathbf{x}}_j$，使用隐变量${\gamma }_{jk},k=1,2...,K$表示来源于哪个成分。
$${\gamma }_{jk}=\{\begin{array}{c}1,{\mathbf{x}}_j来源于第k个成分\\ 0，否则\end{array}$$
完全数据为：$({\mathbf{x}}_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jk}),j=1,2,...m$
令 θ k = { μ k , Σ k , α k } \theta_{k}=\{\bm{\mu}_{k},\bm{\Sigma}_{k},\alpha_{k}\} θk​={μk​,Σk​,αk​}， θ = { θ 1 , θ 2 , . . . , θ K } \bm{\theta}=\{\theta_{1},\theta_{2},...,\theta_{K}\} θ={θ1​,θ2​,...,θK​}
在完全数据情况下的似然函数为：

取对数似然：

E步： 估计隐变量${\gamma }_{jk},,k=1,2,...,K$
求Q函数：

上述式子中的$\log {\alpha }_k$和括号中的内容对于带入的不同的$\gamma$是相同的，所以期望可有上述的代入方式。

M步

下面求解目标函数：
$$\begin{gathered} {\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)}) \\ s.t.\sum _{k=1}^K{\alpha }_k=1 \end{gathered}$$
可以使用拉格朗日乘数法进行求解：
$$L(\theta ,\lambda )=Q(\theta ,{\theta }^{(i)})+\lambda (\sum _{k=1}^K{\alpha }_k-1)$$



类簇划分：
$${\lambda }_j=\arg \underset{k\in 1,2,...,K}{\max }{\hat{\gamma }}_{jk}$$

整体算法的伪代码：