09_期望极大法EM2_统计学习方法

EM内容较多，方便阅读，分成2个部分
上接：09_期望极大法EM1_统计学习方法

四、高斯混合模型GMM

高斯混合模型是指具有如下形式的概率分布模型：
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)\text{\hspace{1em}(20)}$$

• ${\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\left(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\text{\hspace{1em}(21)}$$
  称为第k个分模型。

1、GMM与Adaboost比较

看到式（20）首先让我想到的是Adaboost算法，因为太像了。说下异同点：

• 模型：都是加法模型，Adaboost基学习器不指定，也就是说明也可以使用高斯分布。
• 策略：Adaboost采用的是指数损失函数，而GMM采用的是对数损失函数，根据基学习函数和算法来选择较合适的就行，不是重要的区别项。
• 算法：Adaboost前向分步算法，基学习器数量M是根据终止循环条件确定的，基学习器权重系数${\alpha }_m$之和不指定为1。GMM是采用EM算法，分模型数量K是在学习之前就指定的，分模型权重系数${\alpha }_k$之和指定为1。Adaboost用于监督学习，GMM用于非监督学习，这是两种算法最主要的区别。

2、GMM参数估计的EM算法

假设观测数据$y_1,y_2,\cdots ,y_N$由高斯混合模型生成，
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)\text{\hspace{1em}(22)}$$
其中，$\theta =({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_K;{\theta }_1,{\theta }_2,\cdots ,{\theta }_K)$。

（1）明确隐变量，写出完整数据的对数似然函数

可以设想观测数据$y_j，j=1,2,\cdots ,N$，是这样产生的：首先依概率${\alpha }_k$选择第k个高斯分布分模型$\varphi (y|{\theta }_k)$；然后依第k个分模型的概率分布$\varphi (y|{\theta }_k)$生成观测数据$y_j$。这时观测数据$y_j,j=1,2,\cdots ,N$，是已知的；反映观测数据$y_j$来自第k个分模型的数据是未知的，$k=1,2,\cdots ,K$，以隐变量${\gamma }_{jk}$表示，其定义如下：
$$\begin{gathered} {\gamma }_{jk}=\{\begin{array}{ll}1,& \text{第j个观测来自第k个分模型}\\ 0,& \text{否则}\end{array} \\ j=1,2,\cdots ,N;k=1,2,\cdots ,K \end{gathered}$$
${\gamma }_{jk}$是0-1随机变量。

有了观测数据$y_j$及未观测数据${\gamma }_{jk}$，那么完全数据是
$$(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jK}),j=1,2,\cdots ,N$$
于是，可以写出完全数据的似然函数：
$$\begin{array}{rl}P(y,\gamma |\theta )& =\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},\cdots ,{\gamma }_{jK}|\theta )\\ & =\prod _{j=1}^N\prod _{k=1}^K[{\alpha }_k\varphi (y_j|{\theta }_k){]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\varphi (y_j|{\theta }_k){]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}exp\left(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\end{array}\text{\hspace{1em}(23)}$$
其中，$n_k={\sum }_{j=1}^N{\gamma }_{jk},{\sum }_{k=1}^K{n}_k=N$。

那么，完全数据的对数似然函数为
$$lnP(y,\gamma |\theta )=\sum _{k=1}^K{n}_{k}ln{\alpha }_k+\sum _{k=1}^K\sum _{j=1}^N{\gamma }_{jk}\left[ln\left(\frac{1}{\sqrt{2\pi }}\right)-ln{\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\text{\hspace{1em}(24)}$$

（2）EM算法的E步：确定Q函数

$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E[lnP(y,\gamma |\theta )|y,{\theta }^{(i)}]\\ & =E\left\{\sum _{k=1}^K{n}_{k}ln{\alpha }_k+\sum _{k=1}^K\sum _{j=1}^N{\gamma }_{jk}\left[ln\left(\frac{1}{\sqrt{2\pi }}\right)-ln{\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\right\}\\ & =\sum _{k=1}^K\left\{\sum _{j=1}^N(E{\gamma }_{jk})ln{\alpha }_k+\sum _{j=1}^N(E{\gamma }_{jk})\left[ln\left(\frac{1}{\sqrt{2\pi }}\right)-ln{\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\right\}\end{array}\text{\hspace{1em}(25)}$$
这里需要计算$E({\gamma }_{jk}|y,\theta )$，记为${\hat{\gamma }}_{jk}$。
$$\begin{array}{rl}{\hat{\gamma }}_{jk}& =E({\gamma }_{jk}|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,2,\cdots ,N;k=1,2,\cdots ,K\end{array}\text{\hspace{1em}(26)}$$
${\hat{\gamma }}_{jk}$是在当前模型参数下第$j$个观测数据来自第$k$个分模型的概率，称为分模型$k$对观测数据$y_j$的响应度。

将${\hat{\gamma }}_{jk}=E{\gamma }_{jk}$及$n_k={\sum }_{j=1}^{N}E{\gamma }_{jk}$代入式（25）即得
$$Q(\theta ,{\theta }^{(i)})=\sum _{k=1}^K{n}_{k}ln{\alpha }_k+\sum _{k=1}^K\sum _{j=1}^N{\hat{\gamma }}_{jk}\left[ln\left(\frac{1}{\sqrt{2\pi }}\right)-ln{\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]\text{\hspace{1em}(27)}$$

（3）EM算法的M步

迭代的M步是求函数$Q(\theta ,{\theta }^{(i)})$对$\theta$的极大值，即求新一轮迭代的模型参数：
$${\theta }^{(i+1)}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\text{\hspace{1em}(28)}$$
用${\hat{\mu }}_k,{\hat{\sigma }}_k^2$及${\hat{\alpha }}_k,k=1,2,\cdots ,K$，表示${\theta }^{(i+1)}$的各参数。求${\hat{\mu }}_k,{\hat{\sigma }}_k^2$只需将式（24）分别对${\mu }_k,{\sigma }_k$求偏导数并令其为0，即可得到；求${\hat{\alpha }}_k$是在${\sum }_{k=1}^K{\alpha }_k=1$条件下求偏导数并令其为0得到的。综合，可以构造Q函数拉格朗日函数$L(\theta ,{\alpha }_k)=Q(\theta ,{\theta }^{(i)})+\lambda ({\sum }_{k=1}^K{\alpha }_k-1)$,直接对${\mu }_k,{\sigma }_k，{\alpha }_k$求偏导数并令其为0得到${\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}},k=1,2,\cdots ,K\text{\hspace{1em}(29)}$$
$${\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}(y_j-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K\text{\hspace{1em}(30)}$$
$${\hat{\alpha }}_k=\frac{n_k}{N}=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,\cdots ,K\text{\hspace{1em}(31)}$$
重复以上E步和M步计算，直到对数似然函数值不再有明显的变化为止。

（4）GMM参数估计的EM算法流程

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

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

1）取参数的初始值开始迭代

2）E步：依据当前模型参数，计算分模型k对观测数据$y_j$的响应度
$${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)},j=1,2,\cdots ,N;k=1,2,\cdots ,K$$

3）M步：计算新一轮迭代的模型参数
$${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K$$
$${\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}(y_j-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2,\cdots ,K$$
$${\hat{\alpha }}_k=\frac{n_k}{N}=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2,\cdots ,K$$

4）重复第2）和3）步，直到收敛。

这部分的总结几乎就是照搬李航统计学习方法的内容，刚开始看的时候比较生涩，另外里面的公式（9.28）和（9.29）个人感觉是略有问题，少个求和符号。仔细分析推导后感觉还是挺不错的，所以总结还是统计学习方法内容为主。前面推导EM算法的时候有介绍Andrew NG的推导过程，如果GMM按照这种思路推导也是可以的，可以参照高斯混合模型（GMM），写的非常不错，里面还有详细的的拉格朗日函数求解过程，我也是看过这篇博客后才进一步理解李航的推导过程。

五、EM算法的推广

EM算法还可以解释为F函数的极大极大算法，基于这个解释有若干变形与推广，如广义期望极大算法（GEM）。

1、F函数的极大极大算法

（1）F函数定义

F函数定义：假设隐变量数据Z的概率分布为$\tilde{P}(Z)$，定义分布$\tilde{P}$与参数$\theta$的函数$F(\tilde{P},\theta )$如下：
$$F(\tilde{P},\theta )=E_{\tilde{P}}[lnP(Y,Z|\theta )]+H(\tilde{P})\text{\hspace{1em}(32)}$$
称为F函数。式中$H(\tilde{P})=-E_{\tilde{P}}ln\tilde{P}(Z)=-{\sum }_Z\tilde{P}(Z)ln\tilde{P}(Z)$是分布$\tilde{P}(Z)$的熵。

则F函数可以写作：
$$\begin{array}{rl}F(\tilde{P},\theta )& =E_{\tilde{P}}[lnP(Y,Z|\theta )]+H(\tilde{P})\\ & =\sum _Z\tilde{P}(Z)lnP(Y,Z|\theta )-\sum _Z\tilde{P}(Z)ln\tilde{P}(Z)\\ & =\sum _Z\tilde{P}(Z)ln\frac{P(Y,Z|\theta )}{\tilde{P}(Z)}\\ & \le \sum _{Z}ln\left(\tilde{P}(Z)\frac{P(Y,Z|\theta )}{\tilde{P}(Z)}\right)\\ & =\sum _{Z}lnP(Y,Z|\theta )\\ & =lnP(Y|\theta )\end{array}\text{\hspace{1em}(33)}$$
因为${\sum }_Z\tilde{P}(Z)=1$和$(lnx{)}^{\prime \prime }=-\frac{1}{x^2}\le 0$，所以可以使用Jensen不等式进行上述变形。

之前的EM算法是求期望和求极大值两个步骤。现在的F函数变成了求两次极大值，实际上第一次求极大值也是求期望。把式（10.1）和式（33）的变换过程放在一起看就明了F函数定义，实际上F函数就是对数似然函数的下界。

（2）F函数第1次求极大

F第一次求极大，即使式（33）的变化过程中取得等号，根据Jensen不等式得知，当$\frac{P(Y,Z|\theta )}{\tilde{P}(Z)}=c$,$c$为常数的时候等号成立。根据式（10.2）求解，可得
$$\tilde{P}(Z)=P(Z|Y,\theta )\text{\hspace{1em}(34)}$$

（3）F函数第2次求极大

当求完第1次极大值时可得
$$F(\tilde{P},\theta )=lnP(Y|\theta )$$

所以有
$$\begin{array}{r}{\theta }^{(i+1)}=arg\underset{\theta }{\max }F(\tilde{P},\theta )=arg\underset{\theta }{\max }lnP(Y|\theta )=arg\underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})\end{array}\text{\hspace{1em}(35)}$$

第2步极大值求解和EM算法中的M步一样，个人感觉F函数的极大极大求解和EM就是完全一样，换了个马甲而已。由于第二次极大化求解有时比较困难，所以发展出了不同的GEM算法。

2、GEM算法

（1）GEM算法1

输入：观测数据，F函数；

输出：模型参数.

1）初始化参数${\theta }^{(0)}$,开始迭代

2）第$i+1$次迭代，第一步：记${\theta }^{(i)}$为参数$\theta$的估计值，${\tilde{P}}^{(i)}$为函数$\tilde{P}$的估计，求${\tilde{P}}^{(i+1)}$使$\tilde{P}$极大化$F(\tilde{P},{\theta }^{(i)})$

3）第二步：求${\theta }^{(i+1)}$使$F({\tilde{P}}^{(i+1)},{\theta }^{(i)})$极大化

4）重复2和3，直到收敛

在GEM算法1中，有时求$Q$函数的极大化是很困难的，下面的算法2和算法3并不是求${\theta }^{(i)}$使$Q$函数达到极大的$\theta$，而是找一个${\theta }^{(i+1)}$，使得$Q({\theta }^{(i+1)},{\theta }^{(i)})>Q({\theta }^{(i)},{\theta }^{(i)})$

（2）GEM算法2

输入：观测数据，Q函数；

输出：模型参数.

1）初始化参数${\theta }^{(0)}$,开始迭代

2）第$i+1$次迭代，第一步：记${\theta }^{(i)}$为参数$\theta$的估计值，计算
$$Q(\theta ,{\theta }^{(i)})=E_Z[lnP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _{Z}P(Z|Y,{\theta }^{(i)})lnP(Y,Z|\theta )$$

3)第2步：求${\theta }^{(i+1)}$使
$$Q({\theta }^{(i+1)},{\theta }^{(i)})>Q({\theta }^{(i)},{\theta }^{(i)})$$

4）重复2和3，直到收敛

当参数$\theta$的维度为$d（d\ge 2)$时，可采用一种特殊的GEM算法，它将EM算法的M步分解为d次条件极大化，每次只改变参数向量的一个分量，其余分量不变。

（3）GEM算法3

输入：观测数据，Q函数；

输出：模型参数.

1）初始化参数${\theta }^{(0)}=({\theta }_1^{(0)},{\theta }_2^{(0)},...,{\theta }_d^{(0)})$,开始迭代

2）第$i+1$次迭代，第1步：记${\theta }^{(i)}=({\theta }_1^{(i)},{\theta }_2^{(i)},...,{\theta }_d^{(i)})$为参数$\theta =({\theta }_1,{\theta }_2,...,{\theta }_d)$的估计值，计算
$$Q(\theta ,{\theta }^{(i)})=E_Z[lnP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _{Z}P(Z|Y,{\theta }^{(i)})lnP(Y,Z|\theta )$$

3）第2步：进行d次条件极大化：

首先，在${\theta }_2^{(i)},...,{\theta }_k^{(i)}$保持不变的条件下求使得$Q(\theta ,{\theta }^{(i)})$达到极大的${\theta }_1^{(i+1)}$
然后，在${\theta }_1={\theta }_1^{(i+1)},{\theta }_j={\theta }_j^{(i)},j=3,4,...,k$的条件下求使$Q(\theta ,{\theta }^{(i)})$达到极大的${\theta }_2^{(i+1)}$

如此继续，经过d次条件极大化，得到${\theta }^{(i+1)}=({\theta }_1^{(i+1)},{\theta }_2^{(i+1)},...,{\theta }_d^{(i+1)})$使得
$$Q({\theta }^{(i+1)},{\theta }^{(i)})>Q({\theta }^{(i)},{\theta }^{(i)})$$

4）重复2和3，直到收敛

$\ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast$
未完待续……

参考资料：

• 统计学习方法——李航
• 如何感性地理解EM算法？
• 从最大似然到EM算法浅解
• EM算法学习
• 高斯混合模型（GMM）