从EM算法的典型应用GMM说起，需要知悉的几个点

从EM算法的典型应用GMM说起，需要知悉的几个点

因为毕设一半的工作量是基于GMM的，代码方面也是用c++复现了一遍GMM模型的参数更新公式。大概知道GMM模型的参数更新公式是基于EM算法得出的。但之前一直没有把GMM模型和EM算法的联系真正建立起来，一来也是懒，二来是以前对贝叶斯条件概率那一套没有很好理解。

其实将GMM模型用条件概率的思想进行理解，这就是GMM模型的贝叶斯理解，而利用贝叶斯条件概率得出的生成模型，如果含有隐变量，就自然可以用EM算法来解，（是不是联想到了LDA）（当然有凸函数的限制），而EM算法能保证训练参数的收敛是靠数学上凸函数的Jensen不等式来证明的。

下面大致从6个方面来复习EM算法和复习GMM

• 高斯混合模型（GMM）的表现形式
• GMM的应用
• GMM的贝叶斯理解
• EM算法估计给出的GMM参数更新闭式表达式。
• EM算法收敛性：Jensen不等式
• kmeans和GMM的潜在联系，附带kmeans的优化

高斯混合模型（GMM）的表现形式

设有随机变量$\mathbf{X}$，则混合高斯模型可以用下式表示：
$$p(x)=\sum _{k=1}^K{\pi }_{k}N(x\mid {\mu }_k,{\mathbf{\Sigma }}_k)$$
其中$N(x\mid {\mu }_k,{\mathbf{\Sigma }}_k)$称为混合模型中的第k个分量（component）。
${\pi }_k$是混合系数（mixture coefficient）
满足${\sum }_{k=1}^K{\pi }_k=1$ ,$0<={\pi }_k<=1$ ,
实际上，可以认为${\pi }_k$就是每个分量$N(x\mid {\mu }_k,{\mathbf{\Sigma }}_k)$的权重。

GMM的应用

GMM常用于聚类。如果要从 GMM 的分布中随机地取一个点的话，实际上可以分为两步：首先随机地在这 K 个 Component 之中选一个，每个 Component 被选中的概率实际上就是它的系数${\pi }_{\mathbf{k}}$，选中 Component 之后，再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为已知的问题。

将GMM用于聚类时，假设数据服从混合高斯分布（Mixture Gaussian Distribution），那么只要根据数据推出 GMM 的概率分布来就可以了；然后 GMM 的 K 个 Component 实际上对应K个 cluster 。根据数据来推算概率密度通常被称作 density estimation 。特别地，当我已知（或假定,这里就是假设高斯分布，LDA假设狄利克雷分布）概率密度函数的形式，而要估计其中的参数的过程被称作『参数估计』。

假设有两个聚类，可以定义K=2. 那么对应的GMM形式如下：
$$p(x)={\pi }_{1}N(x\mid {\mu }_1,{\mathbf{\Sigma }}_1)+{\pi }_{2}N(x\mid {\mu }_2,{\mathbf{\Sigma }}_2)$$

上式中未知的参数有六个:$({\pi }_1,{\mu }_1,{\mathbf{\Sigma }}_1;{\pi }_2,{\mu }_2,{\mathbf{\Sigma }}_2)$之前提到GMM聚类时分为两步，第一步是随机地在这K=2分量中选一个，每个分量被选中的概率即为${\pi }_1={\pi }_2=0.5$混合系数，表示每个分量被选中的概率是0.5，即从中抽出一个点，这个点属于第一类的概率和第二类的概率各占一半。但实际应用中事先指定${\pi }_k$的值是很笨的做法,只是第一轮得不得这样初始化而已。当问题一般化后，会出现一个问题：当从图2中的集合随机选取一个点，怎么知道这个点是来$自N(x\mid {\mu }_1,{\mathbf{\Sigma }}_1)$还是$N(x\mid {\mu }_2,{\mathbf{\Sigma }}_2)$呢？换言之怎么根据数据自动确定${\pi }_1，{\pi }_2$的值？这就是GMM参数估计的问题。要解决这个问题，可以使用EM算法。通过EM算法，我们可以迭代计算出GMM中的参数：$({\pi }_k,{\mu }_k,{\mathbf{\Sigma }}_k)$

GMM的贝叶斯理解

上面提到的，“先随机地在这K分量中选一个”，“每个分量被选中的概率” ，这就已经明显感觉到了条件概率的味道了，对，这里我们着重改写GMM的形式，方便我们引出GMM背后的条件概率，后验概率这些东西，从而得到GMM的贝叶斯理解。

GMM的原始形式如下：
$$p(x)=\sum _{k=1}^K{\pi }_{k}N(x\mid {\mu }_k,{\mathbf{\Sigma }}_k)\text{\hspace{1em}(1)}$$
前面提到${\pi }_k$可以看成是第k类被选中的概率。我们引入一个新的K维随机变量$\mathbf{z}$.(贝叶斯公式的味道)。但这个随机变量很特殊，是个二值分布$z_k$(1≤k≤K)只能取0或1两个值。$z_k=1$表示第k类被选中的概率，即：$p(z_k)={\pi }_k$.
如果$z_k=0$表示第k类没有被选中的概率。更数学化一点，$z_k$要满足以下两个条件：

z k ∈ { 0 , 1 } , ∑ K z k = 1 z_k \in \{0,1\},\sum_{K}z_k=1 zk​∈{0,1},K∑​zk​=1

对k=2有两类的例子，则$\mathbf{z}$的维数是2.
如果从第一类中取出一个点，则$\mathbf{z}=(1，0)$,
如果从第二类中取出一个点，则$\mathbf{z}=(0，1)$,
$Z_k$是用来表示第k类被选中的随机变量，只能取0和1两个值（随机变量是要取值的）。
所以变量$Z_k=1$应该是表示“第k类被选中”，而非“第k类被选中的概率”。
但$P(Z_k=1)$表示第k类被选中的概率。
$z_k=1$的概率$P(Z_k=1)$是${\pi }_k$，我们记为$p(z_k)={\pi }_k={\pi }_k^{z_1=1}$
假设$z_k$之间是独立同分布的（iid），我们可以写$\mathbf{z}$的联合概率分布形式，就是连乘：
$$p(\mathbf{z})=p(z_1)p(z_2)...p(z_K)=\prod _{k=1}^K{\pi }_k^{z_k}\text{\hspace{1em}(2)}$$.
因为$z_k$只能取0或1，且$\mathbf{z}$中只能有一个$z_k$为1而其它$(j\ne k)$全为0，所以上式是成立的。(?,这是在解释独立iid吗？解释为什么可以连乘是吗？)
以上其实是我认为 GMM的贝叶斯理解中最难的部分。，因为假设看不懂，后面越糊涂

类确定了（条件），我们可以在类下看数据但高斯分布分布，这个叙述其实可以用条件概率来表示了：$$p(\mathbf{x}|z_k=1)=N(\mathbf{x}\mid {\mu }_k,{\mathbf{\Sigma }}_k)$$
即每一类下数据服从高斯分布，进一步还可以写成：
$$p(\mathbf{x}|\mathbf{z})=N(\mathbf{x}\mid {\mu }_k,{\mathbf{\Sigma }}_k{)}^{z_k}\text{\hspace{1em}(3)}$$

上面式2和式3分别给出了$p(\mathbf{z})$ 和$p(\mathbf{x}|\mathbf{z})$的公式，根据条件概率公式，我们得出观测数据的联合分布：
$$\begin{gathered} p(\mathbf{x})=\sum _{\mathbf{z}}p(\mathbf{z})\ast p(\mathbf{x}|\mathbf{z}) \\ =\sum _{\mathbf{z}}(\prod _{k=1}^K{\pi }_k^{z_k}N(\mathbf{x}\mid {\mu }_k,{\mathbf{\Sigma }}_k{)}^{z_k}) \\ =\sum _{K=1}^K{\pi }_{k}N(\mathbf{x}\mid {\mu }_k,{\mathbf{\Sigma }}_k)\text{\hspace{1em}(4)} \end{gathered}$$
上式从第二个等号到第三个等号其实并不难，第二个等号对z求和，其实就是${\sum }_{K=1}^K$,而里面，只要$i\ne k$,则$z_i=0$,所以z_k=0$的项为1，可忽略，最终得到了第三个等式

可以看到GMM模型的(1)式与(4)式有一样的形式，且(4)式中引入了一个新的变量$\mathbf{z}$，通常称为隐含变量（latent variable）。对于一个要聚类的数据来说，『隐含』的意义是：我们知道数据可以分成两类，但是随机抽取一个数据点，我们不知道这个数据点属于第一类还是第二类，它的归属我们观察不到，因此引入一个隐含变量$\mathbf{z}$来描述这个现象。

再把视角拉回贝叶斯，注意到在贝叶斯的思想下，$p(\mathbf{z})$是先验概率， $p(\mathbf{x}|\mathbf{z})$是似然概率，很自然我们会想到求出后验概率$p(\mathbf{z}|\mathbf{x})$：
$$\begin{gathered} \gamma (z_k)=p(z_k=1|\mathbf{x}) \\ = \\ =\frac{p(z_k=1)p(\mathbf{x}|z_k=1)}{p(\mathbf{x},z_k=1)} \\ = \\ =\frac{p(z_k=1)p(\mathbf{x}|z_k=1)}{{\sum }_{j=1}^{K}p(z_j=1)p(\mathbf{x}|z_j=1)} \\ = \\ =\frac{{\pi }_{k}N(\mathbf{x}\mid {\mu }_k,{\mathbf{\Sigma }}_k)}{{\sum }_{j=1}^K{\pi }_{j}N(\mathbf{x}\mid {\mu }_j,{\mathbf{\Sigma }}_j)}\text{\hspace{1em}(5)} \end{gathered}$$
(上面的第二行，贝叶斯定理，关于这一行的分母，很多人会疑问，应该是$p(\mathbf{x},z_k=1)$还是$p(\mathbf{x})$,按照正常的写法，肯定是$p(\mathbf{x})$，但为了强调$z_k$的取值，有的书会写成$p(\mathbf{x},z_k=1)$，比如李航的小蓝书，这里就约定$p(\mathbf{x})$与$p(\mathbf{x},z_k=1)$是等同的。
第三行全概率公式，最后一个等号，结合三四式.

上式中，我们定义了符号$\gamma (z_k)$来表示第k个分量的后验概率。在贝叶斯的观点下，${\pi }_k$可视为$z_k=1$的先验概率.
上述内容改写了GMM的形式，并引入了隐含变量$\mathbf{z}$和已知$\mathbf{x}$后的的后验概率$\gamma (z_k)$,这样做是为了方便我们使用EM算法来估计GMM的参数。其实EM算法 过程中，我们能看到除了GMM中高斯分布的参数在更新，

隐含变量$\mathbf{z}$的先验概率${\pi }_k$和已知$\mathbf{x}$后的后验概率$\gamma (z_k)$也在不断更新。

EM算法估计给出的GMM参数更新闭式表达式。

注意，式8和式9中，N代表点的数量，$\gamma (z_{nk})$表示点$n({\mathbf{x}}_n)$属于聚类k的后验概率，则，$N_k$代表属于第$k$个聚类的点的数量。$u_k$代表属于第$k$个聚类的所有点的加权平均，每个点的权重是$\gamma (z_{nk})$，即这个点属于聚类k的后验概率大小，跟第k个聚类有关。

EM算法估计GMM参数即最大化(8)，(10)和(12)。需要用到(5)，(8)，(10)和(12)四个公式。我们先指定$\mathbf{\pi }$，$\mathbf{\mu }$和$\mathbf{\Sigma }$的初始值，带入(5)中计算出$\gamma (z_{nk})$，然后再将$\gamma (z_{nk})$带入(8)，(10)和(12)，求得${\mathbf{\pi }}_k$，${\mathbf{\mu }}_k$和${\mathbf{\Sigma }}_k$；接着用求得的${\mathbf{\pi }}_k$，${\mathbf{\mu }}_k$和${\mathbf{\Sigma }}_k$；再带入(5)得到新的$\gamma (z_{nk})$，再将更新后的$\gamma (z_{nk})$带入(8)，(10)和(12)，如此往复，直到算法收敛。

好了，看到这里了，GMM和贝叶斯的关系说清楚了，GMM的参数求解用EM框架给出了流程和闭式表达式了，还有什么呢？有没有看到流程的第五个的“收敛”两个字？

我么难道就不怕这样迭代下去不收敛吗？
所以我们接下来分析下EM 算法的收敛性，准确来说或者具体来说，是含隐变量问题求解过程中的参数更新的收敛性问题。

EM算法收敛性：Jensen不等式

上面谈收敛，首先要弄清楚 ，收敛是对谁说的，即是什么东西收敛，其实这里的收敛是针对参数更新过程中的目标函数而言的。
其实我们在上面讨论时，忽略了谈目标函数的问题，
但其实也提到了目标函数的影子，那就是最大化似然函数。
我们的参数更新过程中，似然函数是会不断变大直至收敛，还是一会儿变大，一会儿变小呢？
从数学的角度分析一下，哎，凸优化知识还有五秒钟抵达战场。。。

还是回忆一下凸函数的图吧，

上面如果t=1/2，就很 直观了
$$f(\frac{1}{2}{x}_1+\frac{1}{2}{x}_2)<=\frac{1}{2}f(x_1)+\frac{1}{2}f(x_2)\text{\hspace{1em}(1)}$$
一般的凸函数定义是
$$f(\lambda x_1+(1-\lambda )x_2)<=\lambda f(x_1)+(1-\lambda )f(x_2)\text{\hspace{1em}(2)}$$

jensen不等式，是对凸函数，把上式扩展到了，无穷多点，即有无穷多个$\lambda$

$$f(\sum _{i=1}^M{\lambda }_i{x}_i)<=\sum _{i=1}^M{\lambda }_{i}f(x_i)\text{\hspace{1em}(3)}$$
$$\sum _{i=1}^M{\lambda }_i=1，{\lambda }_i>0$$
公式(3)被称为 Jensen 不等式，它是式(2)的泛化形式。

因为${\sum }_{i=1}^M{\lambda }_i=1，{\lambda }_i>0$,从概率论的角度 来说，如果把${\lambda }_i$看成是取值为$x_i$的离散变量x的概率分布的话，那么就可以写成期望的形式，如下：
$$f(E[x])<=E[f(x)]\text{\hspace{1em}(4)}$$
$E[.]$表示期望。

同样当函数变为凹函数时，上式的不等号只是变了个方向而已，
$$f(\sum _{i=1}^M{\lambda }_i{x}_i)>=\sum _{i=1}^M{\lambda }_{i}f(x_i)\text{\hspace{1em}(5)}$$
$$f(E[x])>=E[f(x)]\text{\hspace{1em}(6)}$$
为什么要这个，因为在EM算法或者GMM算法的利用Jensen不等式证明收敛性用到的是式子5和式子6这两个不等式。因为EM算法或者GMM算法的利用Jensen证明不等式证明收敛性用到的是Log函数的凹凸性。Log函数是凹函数，二阶导恒小于0.Jensen不等式在Log函数的直观表现如下所示，

下面开始证明之旅。
给定训练样本 { x ( 1 ) , x ( 2 ) , , . . . , , x ( m ) } \{x^{(1)},x^{(2)},,...,,x^{(m)}\} {x(1),x(2),,...,,x(m)},样本间独立，我们想找到每个样本的隐含类别，能使得$p(x,z)$最大。于是$p(x,z)$的最大似然估计是：
$$\begin{gathered} l(\theta )=\sum _{i=1}^{m}logp(x;\theta ) \\ \sum _{i=1}^{m}log\sum _{z}p(x,z;\theta ) \end{gathered}$$

第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求$l(\theta )$一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。

EM是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化$l(\theta )$，我们可以不断地建立$l(\theta )$的下界（E步），然后优化下界（M步）。这句话比较抽象，看下面的

$$\begin{gathered} l(\theta )=\sum _{i=1}logp(x^i;\theta )=\sum _{i=1}log\sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta ) \\ =\sum _{i=1}log\sum _{z^{(i)}}{Q}_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})} \\ \ge \sum _{i=1}\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}\text{\hspace{1em}(1)} \end{gathered}$$

(1)到(2)比较直接，就是分子分母同乘以一个相等的函数。(2)到(3)利用了Jensen不等式，考虑到Log函数是凹函数（二阶导数小于0），而且

这个过程可以看作是对$l(\theta )$求了下界。对于$Q_i$的选择，有多种可能，那种更好的？假设$\theta$已经给定，那么$l(\theta )$的值就决定于$Q_i(z^{(i)})$和$p(x^{(i)},z^{(i)})$了。我们可以通过调整这两个概率使下界J不断上升，$Q_i(z^{(i)})$是可变的，不同的$Q_i(z^{(i)})$和$p(x^{(i)},z^{(i)})$对应了不同的下界J),以逼近$l(\theta )$的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率能够等价于$l(\theta )$了。按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：
$\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}=常数c,Q_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta )}{c}$
c为常数，不依赖$z^{(i)}$，又因为，${\sum }_z{Q}_i(z^{(i)})=1$,所以，
$p(x^{(i)},z^{(i)};\theta )=c\ast Q_i(z^{(i)})$
${\sum }_{z}p(x^{(i)},z^{(i)};\theta )={\sum }_{z}c\ast Q_i(z^{(i)})=c\ast {\sum }_z{Q}_i(z^{(i)})=c$
带入$Q_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta )}{c}$
得到：$$\begin{gathered} Q_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta )}{c} \\ =\frac{p(x^{(i)},z^{(i)};\theta )}{{\sum }_{z}p(x^{(i)},z^{(i)};\theta )} \\ =\frac{p(x^{(i)},z^{(i)};\theta )}{p(x^{(i)},\theta )} \\ =p(z^{(i)}|x^{(i)},\theta ) \end{gathered}$$
至此，我们推出了在固定其他参数$\theta$后，最优的$Q_i(z^{(i)})$（也就是让Jensen不等式等号成立的$Q_i(z^{(i)})$）的计算公式就是后验概率，

解决了$Q_i(z^{(i)})$如何选择的问题。这一步就是E步，建立$l(\theta )$的一个比较好比较紧的下界J。

接下来的M步，就是在给定$Q_i(z^{(i)})$,也就是给定了下界函数J后，调整$\theta$，去极大化下界函数J（从而达到间接极大化$l(\theta )$的目的），（在固定$Q_i(z^{(i)})$后，下界还可以调整的更大，靠调整$\theta$）。

那么究竟怎么确保EM收敛？假定${\theta }^{(t)}$和${\theta }^{(t+1)}$是EM第t次和t+1次迭代后的结果。如果我们证明了$l({\theta }^{(t)})\le l({\theta }^{(t+1)})$，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。

其实综合上面的知识就比较好理解了，
从上面的推导，我一步我们用了最优的$Q_i(z^{(i)})$，也就是后验概率$p(z^{(i)}|x^{(i)},\theta )$，能让Jensen不等式等号成立。

$$l(\theta )=\sum _{i=1}\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{Q_i(z^{(i)})}$$
所以：
$$l({\theta }^{(t+1)})\ge \sum _{i=1}\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};{\theta }^{(t+1)})}{Q_i(z^{(i)})}\text{\hspace{1em}(4)}$$
$$\ge \sum _{i=1}\sum _{z^{(i)}}{Q}_i(z^{(i)})log\frac{p(x^{(i)},z^{(i)};{\theta }^{(t)})}{Q_i(z^{(i)})}\text{\hspace{1em}(5)}$$
$$=l({\theta }^{(t)})\text{\hspace{1em}(6)}$$

上式还是可读，不过是不是有疑问，为什么4式子能大于等于，而6式子就只能等于（中间的5式子是普通的参数估计极大化过程，好理解），因为6式子用了最优的$Q_i(z^{(i)})$，也就是后验概率$p(z^{(i)}|x^{(i)},\theta )$，能让Jensen不等式等号成立。而6式最优的$Q_i(z^{(i)})$不一定在下一次迭代中，同样是最优的，所以需要调整，所以不是等号，而是大于等于。

这样就证明了$l({\theta }^{(t)})$会单调增加。一种收敛方法是$l({\theta }^{(t)})$不再变化，还有一种就是变化幅度很小。

再次解释一下（4）、（5）、（6）。首先（4）对所有的参数都满足，而其等式成立条件只是在固定$\theta$，并调整好Q时成立，而第（4）步只是固定Q，调整$\theta$，不能保证等式一定成立。（4）到（5）就是M步的定义，（5）到（6）是前面E步所保证等式成立条件。也就是说E步会将下界拉到与$l(\theta )$一个特定值（这里$\theta$）一样的高度，而此时发现下界仍然可以上升，因此经过M步后，下界又被拉升，但达不到与$l(\theta )$另外一个特定值一样的高度，之后E步又将下界拉到与这个特定值一样的高度，重复下去，直到最大值。

再啰嗦一句，见上图，在${\theta }^t$处我们又多种$Q_i(z^{(i)})$可选，不同的$Q_i(z^{(i)})$对应不同的下界函数J（图中绿色和蓝色就是两种$Q_i(z^{(i)})$选择下不同的下界函数J），我们会选出最优的$Q_i(z^{(i)})$，让这个$Q_i(z^{(i)})$对应的下界函数正好能在${\theta }^t$处和原目标函数$l({\theta }^t)$相等（其他地方不能保证），这样的$Q_i(z^{(i)})$被我们推导出就是当前的后验概率，然后我们对这个下界函数J，通过参数估计，求极大似然，这里的参数就是$\theta$,于是我们来到了一个新的地方${\theta }^{(t+1)}$.并继续我们的更新之路。

总结

如果将样本看作观察值，潜在类别看作是隐藏变量，那么聚类问题也就是参数估计问题，只不过聚类问题中参数分为隐含类别变量和其他参数，这犹如在x-y坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此什么梯度下降方法就不适用了。但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步估计隐含变量，M步估计其他参数，交替将极值推向最大。EM中还有“硬”指定和“软”指定的概念，“软”指定看似更为合理，但计算量要大，“硬”指定在某些场合如K-means中更为实用（要是保持一个样本点到其他所有中心的概率，就会很麻烦）。

另外，EM的收敛性证明方法确实很牛，能够利用log的凹函数性质，还能够想到利用创造下界，拉平函数下界，优化下界的方法来逐步逼近极大值。而且每一步迭代都能保证是单调的。最重要的是证明的数学公式非常精妙，硬是分子分母都乘以z的概率变成期望来套上Jensen不等式，前人都是怎么想到的。