EM算法、GMM、K-means

   （2019.6.30）这篇博客本来是简单的记录了在读《统计学习方法》——EM算法部分的读后感，但最近因为阅读了《百面机器学习》，让我对EM算法有了更进一步的理解，特对博文进行了更新。
  EM算法是一种迭代算法，用于含有隐含变量的概率模型参数的极大似然估计，或者极大后验概率估计。EM算法的每次迭代由两步组成：E步，求expection；M步，求maximization. 此算法也成为期望极大算法(expection maximization)，EM算法。
  EM算法与初值的选取有关，选择不同的初值可能得到不同的参数估计值。一般用 $Y$代表观测随机变量的数据， $Z$表示隐随机变量的数据。 $Y$和 $Z$连在一起称为完全数据，观测数据 $Y$又称为不完全数据。

1. EM 算法的引入

1.1 EM 算法

   假设给定观测数据$Y$，其概率分布是$P(Y|\theta )$，其中$\theta$是需要估计的模型参数，那么不完全数据$Y$的似然函数是$P(Y|\theta )$，对数似然函数是$L(\theta )=logP(Y|\theta )$;假设$Y$和$Z$的联合概率分布是$P(Y,Z|\theta )$，那么完全数据的对数似然函数是$logP(Y,Z|\theta )$.

• EM算法
    EM算法是通过迭代求$L(\theta )=logP(Y|\theta )$的极大似然估计，每次迭代均分为E步，求期望；M步，求极大化，流程如下：
  输入：观测变量$Y$，隐变量数据$Z$，联合分布$P(Y,Z|\theta )$，条件分布$P(Y|Z,\theta )$;
  输出：模型参数$\theta$.
  (1)选择参数的初值${\theta }^{(0)}$，开始迭代；
  (2)E步：记${\theta }^{(i)}$为第$i$次迭代参数$\theta$的估计值，在第$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)}=argma{x}_{\theta }Q(\theta ,{\theta }^{(i)})$$
  (4)重复第(2)步和第(3)步，直到收敛。
  其中函数$Q(\theta ,{\theta }^{(i)})$是EM算法的核心，成为Q函数。
• Q函数定义：完全数据的对数似然函数$logP(Y,Z|\theta )$关于在给定观测数据$Y$和当前参数${\theta }^{(i)}$下对未观测数据$Z$的条件概率分布$P(Z|Y,{\theta }^{(i)})$的期望称为Q函数，即$$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]$$
  $$=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
• 关于EM算法的几点说明：
  步骤（1）参数的初值可以任意选取，但需注意EM算法对初值选取比较敏感；
  步骤（2）E步求$Q(\theta ,{\theta }^{(i)})$.$Q$函数式中$Z$是未观测数据,$Y$是观测数据.
  步骤（3）M步求 $Q(\theta ,{\theta }^{(i)})$的极大化，得到${\theta }^{(i+1)}$,完成一次迭代，每次迭代使似然函数增大或达到局部最大值；
  步骤（4）给出停止迭代的条件，一般是对比较小的正数${\epsilon }_1,{\epsilon }_2$,若满足$$||{\theta }^{(i+1)}-{\theta }^{(i)}||<{\epsilon }_{1}or||Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})||<{\epsilon }_2$$则停止迭代.

1.2 EM算法的推导

1.2.1 《统计学习方法》中的推导

      通过近似求解观测数据的对数似然函数的极大化问题导出EM算法，从而清楚的看到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 ))-logP(Y|{\theta }^{(i)})$$
根据Jensen不等式:(Jensen不等式是许多不等式推导的基础)$$log\sum _j{\lambda }_j{y}_j\ge \sum _j{\lambda }_{j}log{y}_j$$其中${\lambda }_j\ge 0,{\sum }_j{\lambda }_j=1$. 能够推导出$L(\theta )-L({\theta }^{(i)})$的下界：
$$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)})}-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)})}$$令$$B(\theta ,{\theta }^{(i)})\widehat{=}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)})$是$L(\theta )$的一个下界，且$$L({\theta }^{(i)})=B({\theta }^{(i)},{\theta }^{(i)})$$
因此，任何可以使$B(\theta ,{\theta }^{(i)})$增大的$\theta$，也可以增大$L(\theta )$.为了使$L(\theta )$有尽可能大的增长，选择${\theta }^{(i+1)}$使$B(\theta ,{\theta }^{(i)})$达到极大，即$${\theta }^{(i+1)}=argma{x}_{\theta }B(\theta ,{\theta }^{(i)})$$
省去对参数$\theta$来说是常数的项:$${\theta }^{(i+1)}=argma{x}_{\theta }[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)})}]$$
$$=argmax[\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y|Z,\theta )P(Z|\theta )]$$
$$=argmax[\sum _{Z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )]$$
回忆一下期望表达式：$$E(g(Z))=\sum _{Z_k}P(Z=Z_k)g(Z_k)$$所以上式可化为
$${\theta }^{(i+1)}=argma{x}_{\theta }Q(\theta ,{\theta }^{(i)})$$
对应了EM算法中$Q$函数的定义，上式为等价于EM算法的一次迭代，即求$Q$函数及其极大化。EM算法是不断求解下界的极大化逼近求解对数似然函数极大化算法。

1.2.2 另外一种理解方式

  给出的不完全观测数据 { y 1 , y 2 , ⋯ &ThinSpace; , y N } \lbrace y_1,y_2,\cdots,y_N \rbrace {y1​,y2​,⋯,yN​}，极大化不完全观测数据对数似然函数：
$$\sum _{i}logP(y_i;\theta )=\sum _{i}log\sum _{z^{(i)}}P(y_i,z^{(i)};\theta )$$
直接使用Jensen不等式，且$H(z^{(i)})$是隐变量的一个分布，${\sum }_{z^{(i)}}H(z^{(i)})=1$:
$$\sum _{i}log\sum _{z^{(i)}}P(y_i,z^{(i)};\theta )=\sum _{i}log\sum _{z^{(i)}}H(z^{(i)})\frac{P(y_i,z^{(i)};\theta )}{H(z^{(i)})}\ge \sum _i\sum _{z^{(i)}}H(z^{(i)})log\frac{P(y_i,z^{(i)};\theta )}{H(z^{(i)})}$$
当$\frac{P(y_i,z^{(i)};\theta )}{H(z^{(i)})}=const$(常数)的时候，Jensen不等式取等号，在这里将问题转化为极大化对数似然函数的下界${\sum }_i{\sum }_{z^{(i)}}H(z^{(i)})log\frac{P(y_i,z^{(i)};\theta )}{H(z^{(i)})}$，从而极大化对数似然函数。根据分布$H(z^{(i)})$的概率和是1，以及上述不等式等号成立条件，可以得到$H(z^{(i)})=P(z^{(i)};y_i,\theta )$.
  所以EM算法首先计算出在当前给定的观测数据Y以及当前参数$\theta$下的隐变量概率分布$H(z^{(i)})=P(z^{(i)};y_i,\theta )$，然后极大化对数似然函数下界$$\theta =argma{x}_{\theta }\sum _i\sum _{z^{(i)}}H(z^{(i)})log\frac{P(y_i,z^{(i)};\theta )}{H(z^{(i)})}$$

1.3 EM算法的收敛性

  EM算法提供一种近似计算含有隐变量概率模型的极大似然估计的方法。下面给出EM算法收敛性的两个定理。

定理一

  设$P(Y|\theta )$为观测数据的似然函数，${\theta }^{(i)}(i=1,2,\cdot \cdot \cdot )$为EM算法得到的参数估计序列，$P(Y|{\theta }^{(i)})(i=1,2,\cdot \cdot \cdot )$为对应的对数似然函数序列，则$P(Y|{\theta }^{(i)})$是单调递增的，即$$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$$
  证明：由于$$P(Y|\theta )=\frac{P(Y,Z|\theta )}{P(Z|Y,\theta )}$$
取对数有，$$logP(Y|\theta )=logP(Y,Z|\theta )-logP(Z|Y,\theta )$$
由于$$Q(\theta ,{\theta }^{(i)})=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
设$$H(\theta ,{\theta }^{(i)})=\sum _{Z}logP(Z|Y,\theta )P(Z|Y,{\theta }^{(i)})$$
于是可以将对数似然函数写成$$logP(Y|\theta )=Q(\theta ,{\theta }^{(i)})-H(\theta ,{\theta }^{(i)})$$
分别取$\theta$为${\theta }^{(i)}$和${\theta }^{(i+1)}$并相减，有$$logP(Y|{\theta }^{(i+1)})-logP(Y|{\theta }^{(i)})$$ $$=[Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i)},{\theta }^{(i)})]-[H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i)},{\theta }^{(i)})]$$
式子右端的第1项，由于${\theta }^{(i+1)}$使$Q(\theta ,{\theta }^{(i)})$达到极大，所以有$$Q({\theta }^{(i+1)},{\theta }^{(i)})\ge 0$$
右端第二项：
$$H({\theta }^{(i+1)},{\theta }^{(i)})=\sum _Z[log\frac{P(Y|Z,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}]P(Z|y,{\theta }^{(i)})$$
$$\le log[\sum _Z\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}P(Z|Y,{\theta }^{(i)}]$$
$$=log(\sum _{Z}P(Z|Y,{\theta }^{(i+1)}))=0$$定理一得证。

定理二

  设$L(\theta )=logP(Y|\theta )$为观测数据的对数似然函数，${\theta }^{(i)}(i=1,2,\cdot \cdot \cdot )$为EM算法得到的参数估计序列，$L({\theta }^{(i)})(i=1,2,\cdot \cdot \cdot )$为对应的对数似然函数序列.
(1)如果$P(Y|\theta )$有上界，则$L({\theta }^{(i)})=logP(Y|{\theta }^{(i)})$收敛到某一值$L^{\ast }$;
(2)在函数$Q(\theta ,{\theta }^{{}^{\prime }})$与$L(\theta )$满足一定条件下，由EM算法得到的参数估计序列${\theta }^{(i)}$的收敛值${\theta }^{\ast }$是$L(\theta )$的稳定点.

2. 高斯混合模型

2.1 高斯混合模型

• 定义9.2(高斯混合模型) 高斯混合模型是指具有如下形式的概率分布模型：$$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^2)$,$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt[]{2\pi }{\sigma }_k}exp(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2})$$称为第$k$个模型。

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

    假设观测数据$y_1,y_2,\cdots ,y_N$由高斯混合模型生成$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$其中，$\theta =({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K;{\theta }_1,{\theta }_2,\cdots ,{\theta }_K)，{\theta }_i=({\mu }_i,{\sigma }_i^2),\theta$就是高斯混合模型中需要确定的参数。
     对于多个高斯分布混合模型的EM算法理解，可以在单个高斯模型的基础上进行等价考虑。

• 高斯混合模型的EM算法：
      采用EM算法对高斯混合模型参数进行迭代求解，根据当前模型参数，定义分模型$k$对观测数据$y_j$的影响度（当前的$y_j$的百分之多少是来自模型$j$）
  $${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)}$$表示当前模型下第$j$个数据来自第$k$个分模型的概率。
      然后参数${\hat{\mu }}_k、{\hat{\sigma }}_k^2、{\hat{\alpha }}_k$的公式为：$${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}$$
  $${\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}(y_j-{\hat{\mu }}_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}$$
  $${\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N}$$
  重复上面的步骤，直到达到迭代停止条件。

3. K-means

书中例题

   将可观测数据表示为$Y=(Y_1,Y_2,\cdot \cdot \cdot ,Y_n{)}^T$,未观测数据表示为$Z=(Z_1,Z_2,\cdot \cdot \cdot ,Z_n{)}^T$,则观测数据的似然函数为$$P(Y|\theta )=\sum _{Z}P(Z|\theta )P(Y|Z,\theta )$$即$$P(Y|\theta )=\prod _{j=1}^n[\pi p^{yi}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i}]$$
考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，即$$\hat{\theta }=argma{x}_{\theta }logP(Y|\theta )$$这个问题没有解析解，只有通过迭代的方法求解，EM就是可以求解这类问题的一种算法。

• 书中示例的 EM算法求解过程：
    首先选取参数的初值，记作${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$，通过下面步骤迭代计算参数的估计值直到收敛为止。记第$i$次迭代参数的估计值为${\theta }^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$.EM算法的第$i+1$次迭代如下：
   $E$步：计算在模型参数${\pi }^{(i)},p^{(i)},q^{(i)}$下观测数据$y_i$来自掷硬币$B$的概率$${\mu }_j^{(i+1)}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_j}(1-p^{(i)}{)}^{1-y_j}}{{\pi }^{(i)}(p^i{)}^{y_j}(1-p^{(i)}{)}^{1-y_i}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_j}(1-q^{(i)}{)}^{1-y_i}}$$
   $M$步：计算模型参数的新估计值：$${\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i+1)}$$ $$p^{(i+1)}=\frac{{\sum }_{j=1}^n{\mu }_j^{(i+1)}{y}_i}{{\sum }_{j=1}^n{\mu }_j^{(i+1)}}$$ $$q^{(i+1)}=\frac{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})y_j}{{\sum }_{j=1}^n(1-{\mu }_j^{(i+1)})}$$