统计学习方法9—EM算法

---

title: EM算法
date: 2019-07-25 14:12:12
tags: [机器学习,EM算法]
categories: 机器学习

  EM算法是一种迭代算法，是一种用于计算包含隐变量概率模型的最大似然估计方法，或极大后验概率。EM即expectation maximization，期望最大化算法。

1. 极大似然估计

  在概率模型中，若已知事件服从的分布或者其他概率模型的参数，那么我们可以通过计算得到某事件发生的概率。而在估计中，这些变成了方向过程：已知一组数据发生的结果，相当于获得了经验概率，通过这组数据假设模型服从什么分布，再通过经验概率求解模型参数。
  比如统计学校学生身高服从的概率分布，抽样1000人得到他们的身高数据，再假设身高服从正态分布，于是只需要求解$u,\sigma$即可。剩下的问题就是如何确定这两个参数，极大似然估计的思想就是求解在这两个参数下，得到的 这组样本发生的概率最大的情况便是这两个参数最有可能的取值。而抽取每个人的身高都是独立的，让每个p最大化就是让其乘积最大化，于是得到最大似然函数
$$L(\theta )=\prod _{x}p(x|\theta )$$
  再求解其关于$\theta$的极大化问题便能得到一个对$\theta$的估计
$$\hat{\theta }=argmaxL(\theta )$$

2.EM算法

  上面是模型变量都是可观测的情况，这种情况可以直接根据观测到的样本数据，使用极大似然估计对模型参数进行估计。但若模型中含有隐变量的时候，就不能简单的使用极大似然估计进行计算。EM算法就是针对含有隐变量的概率模型参数的极大似然估计方法。
例子：
  假设三枚硬币，记为A、B、C，它们正面向上的概率分别为n、p、q。实验如下：先抛A，根据A的结果选择B或者C再次抛，将这次正面向上的结果记为1，反面记为0，其中A正面向上则选择B，反面则选择C。经过n次实验后得到n个观测值，根据这n个观测值估计n、p、q的参数值。实验过程中只能观测到最后结果，并不能观察实验过程，也就是A的结果是未知的。该模型可以表示如下
$$P(y|\theta )=\sum _{z}P(y,z|\theta )=\sum _{z}P(z|\theta )P(y|z,\theta )$$
  其中y表示随机变量Y的观测值，表示该次结果是1或0。z代表隐变量，表示模型中未能观测到的A的结果，$\theta =（n,p,q)$是模型参数。其中Y是不完全数据，Y+Z是完全数据。
  若是没有隐变量Z，那么可直接使用对数极大似然函数估计
$$\hat{\theta }=arg\underset{\theta }{\max }L(\theta )=arg\underset{\theta }{\max }\sum _{y}logP(y|\theta )$$
  加入隐变量后，对数极大似然函数变为
$$L(\theta ,z)=\sum _{y}log\sum _{z}P(y,z|\theta )$$
  求解
$$\hat{\theta },z=arg\underset{\theta ,z}{\max }L(\theta ,z)$$
  按照极大似然估计中的方法，似乎应该对$\theta ,z$求导然后令其为0解方程组，然后注意此时的$L(\theta ,z)$函数，log内含有求和运算，若是直接求导计算十分困难，因此退而求其次，既然要极大化$L(\theta ,z)$，就先找到一个它的下界，再想办法提高这个下界，同样能达到极大化L的效果。使用Jense不等式对似然函数进行放缩。

• Jense不等式（琴生不等式）

  凸函数：是定义在实数域上的函数，如果对于任意的实数，都有：
$$f^{\prime \prime }⩾0$$
   则该函数是凸函数。
  当函数是凸函数的时候，Jense不等式的含义是函数的期望大于等于期望的函数（凹函数相反）。图下图所示

  二维情况下可用凸函数定义来解释，当一个函数是凸函数的时候，它满足
$$f(\frac{a+b}{2})⩽\frac{f(a)+f(b)}{2}$$
  左边其实相当于其变量x先求期望后带入函数，右边是先求函数值再计算函数值的期望，也就是
$$f(E(x))⩽E(f(x))$$
  再回到$L(\theta ,z)$中来，目的是为了将对数函数中的和项去掉，便可利用$jense$不等式的性质，将期望的函数变为函数的期望。先进行放缩
$$\begin{array}{rl}L(\theta ,z)& =\sum _{y}log\sum _{z}P(y,z|\theta )\\ & =\sum _{y}log\sum _z{Q}_i(z)\frac{P(y,z|\theta )}{Q_i(z)}\\ & ⩾\sum _y\sum _z{Q}_i(z)log\frac{P(y,z|\theta )}{Q_i(z)}\end{array}$$
  其中最后一步用到了$jense$不等式，因为对数函数是凹函数，所以不等号反了过来，$f(E(x))⩾E(f(x))$,此处
$$f(x)=log(x),x=\frac{P(y,z|\theta )}{Q_i(z)}，$$
  并且所添加的$Q_i(x)$满足$$\sum _z{Q}_i(z)=1$$  这是根据第三类$jense$不等式的条件设定的，不同系数的加权求和期望只要满足系数之和为1就能使用$jense$不等式。
  所以得到结论，$log\frac{P(y,z|\theta )}{Q_i(z)}$的加权平均就是$L(\theta ,z)$的一个下界。这便是EM算法中E（期望）的来由。
  目前$Q(z)$还是未知的，需要根据一些条件来选择一个合适的函数，再次强调最终目的是极大化似然函数，现在我们得到了似然函数的一个下界，一个想法就是让这个下界去更好的逼近似然函数的真实值，下界是根据不等式放缩后得到的，那么当不等式取等号的时候便是最接近原函数的时候。回忆$jense$不等式$f(E(x))⩽E(f(x))$，显然当$x$为常数的时候不等式成立。即
$$x=\frac{P(y,z|\theta )}{Q_i(z)}=c$$
  既然要确定的是$Q_i(z)$，不妨设$\theta$为常数，事实上，EM是一种迭代算法，也就是我们是先给出$\theta$的假设值后再进行迭代计算的，设上一次为第i次迭代，${\theta }^{(i)}$为第i次的参数值，由上式得
$$\begin{array}{rl}P(y,z|{\theta }^{(i)})& =c{Q}_i(z)\\ \sum _{z}P(y,z|{\theta }^{(i)})& =\sum _{z}c{Q}_i(z)\\ \sum _{z}P(y,z|{\theta }^{(i)})& =c\end{array}$$
  将c带入x
$$\begin{array}{rl}{Q}_i(z)& =\frac{P(y,z|{\theta }^{(i)})}{{\sum }_{z}P(y,z|{\theta }^{(i)})}\\ & =\frac{P(y,z|{\theta }^{(i)})}{P(y|{\theta }^{(i)})}\\ & =P(z|y,{\theta }^{(i)})\end{array}$$
  第二步用到边缘概率，第三步条件概率。至此，我们得出了Q(z)的表达式，Q(z)是已知样本观测和模型参数下的隐变量的分布。将Q_i(z)视作常数去掉后，得到下面表达式

$$\begin{array}{rl}L(\theta ,z)& =\sum _y\sum _z{Q}_i(z)log\frac{P(y,z|\theta )}{Q(z)}\\ & =\sum _y\sum _z{Q}_i(z)logP(y,z|\theta )\\ & =\sum _y\sum _{z}P(z|y,{\theta }^{(i)})logP(y,z|\theta )\end{array}$$
  到这儿已经完成了E步，对对数似然函数下界的期望的表示，这个期望是$logP(y,z|\theta )$关于$Q_i(z)=P(z|y,{\theta }^{(i)})$的期望，接下来需要做的便是极大化这个期望，也就是M步。极大化的方法就是我们常见的优化问题，使用拉格朗日乘数法即可，添加约束条件${\sum }_{z}P(z|y,{\theta }^{(i)})=1$后对$\theta ,z,\lambda$求偏导令为0。求解

$$\hat{\theta }=arg\underset{\theta }{\max }L(\theta ,z)$$

  迭代的停止条件常常为$||{\theta }^{(i+1)}-{\theta }^{(i)}||⩽\epsilon ,or,||L({\theta }^{(i+1)})-L({\theta }^{(i)}||⩽\epsilon$

3.EM算法的收敛性

  观测数据的似然函数是$L(\theta )=P(y|\theta )$，我们先证明似然函数关于每次迭代后的$\theta$是单调递增的
proof：
$$P(Y|\theta )=\frac{P(Y,Z|\theta )}{P(Z|Y,\theta )}$$
取对数
$$logP(Y|\theta )=logP(Y,Z|\theta )-logP(Z|Y,\theta )$$
由前面结论
$$L(\theta ,z)=\sum _{z}P(Z|Y,{\theta }^{(i)})logP(Y,Z|\theta )$$
令
$$H(\theta ,z)=\sum _{z}P(Z|Y,{\theta }^{(i)})logP(Z|Y,\theta )$$
于是，对数似然函数可写作
$$logP(Y|\theta )=L(\theta ,z)-H(\theta ,z)$$
那么
$$logP(Y|{\theta }^{(i+1)})-logP(Y|{\theta }^{(i)}=[L({\theta }^{(i+1)},z)-L({\theta }^{(i)},z)]-[H({\theta }^{(i+1)},z)-H({\theta }^{(i)},z)]$$
只需证右边非负，因${\theta }^{(i+1)}$使原似然函数$L({\theta }^{(i+1)},z)$达到极大，所以
$$L({\theta }^{(i+1)},z)-L({\theta }^{(i)},z)⩾0$$
对于第二项
$$\begin{array}{rl}& H({\theta }^{(i+1)},z)-H({\theta }^{(i)},z)\\ & =\sum _{z}P(Z|y,{\theta }^{(i)})log\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}\\ (Jensen不等式)& ⩽log\left(\sum _{z}P(Z|Y,{\theta }^{(i)})\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}\right)\\ & =log\sum _{z}P(Z|Y,{\theta }^{(i+1)})\\ & =log1=0\end{array}$$
  因此$logP(Y|\theta )$关于迭代的$\theta$单调递增，那么如果$P(Y|\theta )$有上界，经过$\theta$的单增迭代一定会收敛到某一值，还有一个定理是EM算法得到的$\theta$收敛值是似然函数的稳定点，人生苦短，证明就免了吧。

4.EM算法在高斯混合模型中的参数估计

  通过前面的介绍已经知道EM算法用于解决含隐变量的概率模型，EM算法可以认为是一种解决问题的思想，在不知道隐变量发生情况的条件下，我们通过计算已知观测数据，求出不同数据在隐变量的某个事件下发生的期望，再通过这个期望去最大化该期望下的似然函数，如此迭代下去便能得到一个局部极大的结果。
  高斯混合模型也是这样的一个概率模型，它由多个正态分布加权组合而成。也就是可以观测到混合模型产生的样本的分布，但是不知道这每一个样本是由其中哪个高斯模型产生的。因此，隐变量就是哪个样本选择了哪个模型。

• 高斯混合模型

  假设观测数据$y_1,y_2,...,y_n$由高斯混合模型产生
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\Phi (y|{\theta }_k)$$
  其中$\theta =({\alpha }_1,...,{\alpha }_k;{\theta }_1,...,{\theta }_k)$，由前分析得隐变量是每个样本来自哪个模型，j个样本和k个模型，隐变量${\gamma }_{jk}$表示如下
$${\gamma }_{jk}＝\{\begin{array}{rlr}& 1,& 第j个观测来自第k个分模型\\ & 0,& 否则\end{array}$$
  于是有$j$个观测数据$j\ast k$个位观测数据，总共$(j+1)\ast k$个数据便组成了完全数据，完全数据的似然函数为
$$\begin{array}{rl}P(y,\gamma |\theta )& =\prod _j^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jk}|\theta )\\ & =\prod _k\prod _j[{\alpha }_k\Phi (y_j|{\theta }_k){]}^{{\gamma }_{jk}}\\ & =\prod _k{\alpha }_k^{n_k}\prod _j[\Phi (y_j|{\theta }_k){]}^{{\gamma }_{jk}}\end{array}$$
  其中$n_k={\sum }_{j=1}^N{\gamma }_{jk},{\sum }_{k=1}^K{n}_k=N$,
$$\Phi (y_j|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\delta }_k}exp\left(\frac{(y_i-u_k{)}^2}{2{{\delta }_k}^2}\right)$$
  显然$n_k$表示的是第k个模型出现的样本数量，因此所有模型出现的样本数量和就是N，完全数据的对数似然函数为
$$logP(y,\gamma |\theta )=\sum _k\left\{n_{k}log{\alpha }_k+\sum _j{\gamma }_{jk}[log(\frac{1}{\sqrt{2\pi }}-log{\delta }_k-\frac{(y_j-u_k{)}^2}{2{{\delta }_k}^2})]\right\}$$
  确定好完全数据的对数似然函数后就可以直接套用EM算法的框架了，首先是E-step，求解对数似然函数关于隐变量的后验概率的期望
$$\begin{array}{rl}L(\theta ,\gamma )& =E[logP(y,\gamma |\theta );y,{\theta }^{(i)}]\\ & =E\left\{\sum _k\left\{n_{k}log{\alpha }_k+\sum _j{\gamma }_{jk}[log(\frac{1}{\sqrt{2\pi }}-log{\delta }_k-\frac{(y_j-u_k{)}^2}{2{{\delta }_k}^2})]\right\}\right\}\\ & =\sum _k\left\{\sum _{j}E({\gamma }_{jk})log{\alpha }_k+\sum _{j}E{\gamma }_{jk}[log(\frac{1}{\sqrt{2\pi }}-log{\delta }_k-\frac{(y_j-u_k{)}^2}{2{{\delta }_k}^2})]\right\}\end{array}$$
  注意其中的$n_k$是关于${\gamma }_{jk}$的函数$n_k={\sum }_j{\gamma }_{jk}$，注意上面两步的化简中，并没有直接计算期望，而是根据期望的性质对式子进行化简，将除了$\gamma$的式子看作常数，由$E(aX)=aE(x)$等期望的运算法则化简到只需要求$E({\gamma }_{jk})$。
  计算$E({\gamma }_{jk})$,记为$\hat{\gamma }$，注意这个期望是对$P({\gamma }_{jk}|y,\theta )$的期望。
$$\begin{array}{rl}\hat{\gamma }& =E({\gamma }_{jk}|y,\theta )\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )}{{\sum }_{k}P({\gamma }_{jk}=1,y_j|\theta )}\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )P({\gamma }_{jk}=1|\theta )}{{\sum }_{k}P({\gamma }_{jk}=1,y_j|\theta )P({\gamma }_{jk}=1|\theta )}\\ & =\frac{{\alpha }_k\Phi (y_j|{\theta }_k)}{{\sum }_k{\alpha }_k\Phi (y_j|{\theta }_k)}\end{array}$$
  带入$L(\theta ,\gamma )$得
$$L(\theta ,\gamma )=\sum _k\left\{n_{k}log{\alpha }_k+\sum _j{\hat{\gamma }}_{jk}[log(\frac{1}{\sqrt{2\pi }}-log{\delta }_k-\frac{(y_j-u_k{)}^2}{2{{\delta }_k}^2})]\right\}$$
  接下就是M步，只需最大化这个L即可
$${\theta }^{(i+1)}=arg\underset{\theta }{\max }L(\theta ,\gamma )$$
  求解方式同样是拉格朗日，约束条件${\sum }_k{\alpha }_k=1$