《统计学习方法》第九章: EM算法及其推广 读书笔记

9.EM算法及其推广

• EM算法(expectation maximization algorithm,期望极大算法)
• 是一种非监督模型
• 是含有缺失数据的概率模型参数的极大似然估计法
• 算法每次迭代分两步：
  • E：求期望
  • M：求极大

其实k-means聚类和高斯混合模型都是EM的推广。
分类模型试图从数据的内在联系分析出数据可以分为几类，分别属于哪一类。

9.1概念

• 不完全数据：观测数据$X$(观测随机变量得到的结果数据)
• 完全数据：观测数据$X$和隐随机变量的数据 $Z$。
  我的理解，对于分类模型，这里的隐随机变量就是数据的类别。

9.1.1模型解释

设待估计的模型参数为$\theta$。
例如对于k-means来说，$\theta$就是各聚类的中心${\mu }_1,\cdots ,{\mu }_k$；隐变量Z就是最终的K个分类$1,\cdots ,k$。
对于混合高斯分布来说，$\theta$就是各高斯分布的参数${\alpha }_i,{\mu }_i,{\Sigma }_i$；隐变量Z就是K个分布$1,\cdots ,k$。

• 每个样本$x_i$的真实类别$z_i$是隐随机变量，未知；所以EM算法的步骤：
  • 初始化${\theta }^0$
  • E步: 计算$E(z_i)$，$E(z_i)$可用${\theta }^{(n)}$表示。即计算在参数值为${\theta }^{(n)}$的情况下，样本真实类别的期望$E(z_i)$。
    对于k-means，这一步计算的是在当前聚类中心为${\mu }_1^{(n)},\cdots ,{\mu }_k^{(n)}$的条件下，样本的可能分类${\hat{z}}_i$。
  • M步：用$E(z_i)$代替$z_i$带入$L(\theta )$，求本轮迭代中使得极大似然函数最大的$\theta$，即
    ${\theta }^{(n+1)}=\arg {\max }_{\theta }L(\theta )$。
    对于k-means来说，即按照上一轮聚类中心将样本集划分后，将聚类中心更新，值为当前分类子集的质心。

9.1.2极大似然函数

当数据完整时

• $X$和$Z$的联合概率分布为$P(x,z|\theta )$
• 极大似然函数
  $$P(X,Z|\theta )=\prod _{i=1}^{m}P(x_i,z_i|\theta )$$
• 对数极大似然函数为
  $$L(\theta )=\log P(X,Z|\theta )=\log \prod _{i=1}^{m}P(x_i,z_i|\theta )=\sum _{i=1}^m\log P(x_i,z_i|\theta )$$

当数据不完整时

• 极大似然函数，假设数据集共有m个样本
  $$l=P(X|\theta )=\prod _{i=1}^{m}P(x_i|\theta )=\prod _{i=1}^m(\sum _{z_i=j}^{k}P(x_i,z_i|\theta ))$$
  其中$j=1,\cdots ,k$
• 对数极大似然函数
  $$L(\theta )=\log P(X|\theta )=\sum _{i=1}^m\log P(x_i|\theta )=\sum _{i=1}^m\log \sum _{z_i=j}^{k}P(x_i,z_i|\theta )$$

9.2EM算法

9.2.1推导

EM算法通过一步步迭代$\theta$值，逐步最大化似然函数。
假设在第n次迭代取值${\theta }^{(n)}$，更新$\theta$值时，希望$L(\theta )-L({\theta }^{(n)})>0$，以此逐步最大化似然函数。
考虑两者的差
$$\begin{array}{rl}L(\theta )-L({\theta }^{(n)})& =\sum _{i=1}^m\log \sum _{z_i=j}^{k}P(x_i,z_i|\theta )-\sum _{i=1}^m\log P(x_i|{\theta }^{(n)})\\ & =\sum _{i=1}^m\log [\sum _{z_i=j}^k(P(z_i|x_i,{\theta }^{(n)})\frac{P(x_i,z_i|\theta )}{P(z_i|x_i,{\theta }^{(n)})})]-\sum _{i=1}^m\log P(x_i|{\theta }^{(n)})\end{array}$$

• 有Jensen不等式，当$f(x)$是凸函数时$f({\sum }_i{\alpha }_i{x}_i)⩾{\sum }_i{\alpha }_{i}f(x_i)$；此时$f(x)=\log x$，满足不等式
• 且${\sum }_{z_i=j}^{k}P(z_i|x)=1$

根据两个性质，得到：
$$\begin{array}{rl}L(\theta )-L({\theta }^{(n)})& =\sum _{i=1}^m\log [\sum _{z_i=j}^k(P(z_i|x_i,{\theta }^{(n)})\frac{P(x_i,z_i|\theta )}{P(z_i|x_i,{\theta }^{(n)})})]-\sum _{i=1}^m\log P(x_i|{\theta }^{(n)})\\ & ⩾\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log \frac{P(x_i,z_i|\theta )}{P(z_i|x_i,{\theta }^{(n)})}-\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log P(x_i|{\theta }^{(n)})]\\ & =\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log \frac{P(x_i,z_i|\theta )}{P(z_i|x_i,{\theta }^{(n)})P(x_i|{\theta }^{(n)})}]\\ & =\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log \frac{P(x_i,z_i|\theta )}{P(x_i,z_i|,{\theta }^{(n)})}]\\ \text{令}B(\theta ,{\theta }^{(n)})& =L({\theta }^{(n)})+\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log \frac{P(x_i,z_i|\theta )}{P(x_i,z_i|,{\theta }^{(n)})}]\\ \text{则}L(\theta )& ⩾B(\theta ,{\theta }^{(n)})\end{array}$$
上式当$\theta$取${\theta }^{(n)}$时等号成立,证：
$$\begin{array}{rl}B({\theta }^{(n)},{\theta }^{(n)})& =L({\theta }^{(n)})+\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log \frac{P(x_i,z_i|{\theta }^{(n)})}{P(x_i,z_i|,{\theta }^{(n)})}]\\ & =L({\theta }^{(n)})\end{array}$$
此时$B(\theta ,{\theta }^{(n)})$相当于$L(\theta )$的下界，如果能最大化$B(\theta ,{\theta }^{(n)})$，也能够使$L(\theta )$增大。
更新${\theta }^{(i+1)}$为使$B(\theta ,{\theta }^{(i)})$最大的值
$$\begin{array}{rl}{\theta }^{(n+1)}& =\arg \underset{\theta }{\max }B(\theta ,{\theta }^{(n)})\\ & =\arg \underset{\theta }{\max }L({\theta }^{(n)})+\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log \frac{P(x_i,z_i|\theta )}{P(x_i,z_i|,{\theta }^{(n)})}]\\ & =\arg \underset{\theta }{\max }\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log P(x_i,z_i|\theta )]+constant\\ & =\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(n)})\end{array}$$
其中$Q(\theta ,{\theta }^{(n)})={\sum }_{i=1}^m[{\sum }_{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log P(x_i,z_i|\theta )]={\sum }_{i=1}^{m}E(\log P(x_i,z_i|\theta )|x_i,{\theta }^{(n)})$;
其中$P(z_i|x_i,{\theta }^{(n)})$是在给定观测数据$x_i$和当前参数${\theta }^{(n)}$下，对隐变量$z_i$的期望。
下图为不完全数据的对数似然函数$L(\theta )$和$B(\theta ,{\theta }^{(i)})$的关系

当$L(\theta )$取${\theta }^{(i)}$时，求$B(\theta ,{\theta }^{(i)})$曲线，求得时$B(\theta ,{\theta }^{(i)})$最大点${\theta }^{(i+1)}$；令$L(\theta )$取${\theta }^{(i+1)}$时，继续求$B(\theta ,{\theta }^{(i+1)})$曲线，继续循环。

9.2.2算法

(EM算法选择不同的初值可能得到不同的参数估计值)

输入：观测变量X，隐变量数据Z，联合分布$P(X,Z|\theta )$,条件分布$P(Z|X,\theta )$(求得Z的期望，带入联合分布中，以求未知数的最大值)
输出：模型参数$\theta$

• 1)选择参数的初始值${\theta }^{(0)}$开始迭代；
• 2)E步：${\theta }^{(n)}$为第n次迭代参数的估计值，在第n+1次迭代的E步，计算
  $$\begin{array}{r}Q(\theta ,{\theta }^{(n)})=\sum _{i=1}^m[\sum _{z_i=j}^{k}P(z_i|x_i,{\theta }^{(n)})\log P(x_i,z_i|\theta )]\end{array}$$
  此时$P(z_i|x_i,{\theta }^{(n)})$就是在给定观测数据$x_i$和当前参数估计${\theta }^{(n)}$下隐变量$z_i$的条件概率分布
• 3)M步:求使$Q(\theta ,{\theta }^{(n)})$最大化的$\theta$，确定第n+1次的参数估计值${\theta }^{(n+1)}$
  $${\theta }^{(n+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(n)})$$
• 4)重复2、3步，直到收敛(停止的条件，例如$||{\theta }^{(n+1)}-{\theta }^{(n)}||⩽{ϵ}_1$)

9.2.3 算法的收敛性

• 观测数据的似然函数$P(X|\theta )$在$\theta$的迭代过程中单调递增
• EM算法得到的$\theta$值是观测数据对数似然函数$L(\theta )=\log P(X|\theta )$的局部最优解

9.3EM在高斯混合模型中的应用

高斯混合模型：样本以不同的可能性来自不同的高斯分布。
$$\begin{array}{rl}P(x|\theta )& =\sum _{z}p(x,z|\theta )=\sum _{z}p(z|\theta )p(x|z,\theta )\\ & =\sum _{k=1}^K{\alpha }_k\varphi (x|{\theta }_k)\end{array}$$

• 其中$z=1,\cdots ,K$是隐数据，代表取自第几个高斯分布，其中$p(z=k|\theta )={\alpha }_k$，代表第k个分模型的权重，${\alpha }_k⩾0,{\sum }_{k=1}^K{\alpha }_k=1$；

• $p(x|z=k,\theta )=\varphi (x|{\theta }_k)$，$\varphi (x|{\theta }_k)$是高斯分布密度，${\theta }_k=({\mu }_k,{\sigma }_k^2)$,其中第k个分模型
  $$\varphi (x|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}exp(-\frac{(x-{\mu }_k{)}^2}{2{\sigma }_k^2})$$

• $\theta =({\alpha }_1,\cdots ,{\alpha }_K;{\theta }_1,\cdots ,{\theta }_K,{\sigma }_1^2,\cdots ,{\sigma }_K^2,)$

• 观测数据x，$x_j=1,2,\cdots ,N$

• 隐变量${\gamma }_{jk}$
  $${\gamma }_{jk}=\{\begin{array}{ll}1,& \text{第j个观测来自第k个分模型}\\ 0,& \text{否则}\end{array}$$
  因此$E{\gamma }_{jk}=P({\gamma }_{jk}=1)$
  (样本中只有观测数据$x_j$，并不知道其由哪个模型生成的观测数据)

• 观测数据$x_j$对应的未观测数据$r_j=(r_{j1},r_{j2},\cdots ,r_{jK})$，取值可能为$(1,0,\cdots ,0),(0,1,\cdots ,0),\cdots ,(0,0,\cdots ,1)$,只有一个值为1，其他都为0

9.3.1利用EM结论的推导

EM算法的Q函数
$$\begin{array}{r}Q(\theta ,{\theta }^{(n)})=\sum _{i=1}^m[\sum _{z_i=j}^{K}P(z_i|x_i,{\theta }^{(n)})\log P(x_i,z_i|\theta )]\end{array}$$
对于高斯混合模型来说:(其中${\theta }_{z_i}$是总参数$\theta$中第$z_i$个高斯分布的参数，${\alpha }_{z_i}$是第$z_i$个高斯分布的权重)
$$\begin{array}{rl}P(x_i,z_i|\theta )& =P(x_i|z_i,\theta )P(z_i|\theta )={\alpha }_{z_i}\varphi (x_i|{\theta }_{z_i})\\ P(z_i|x_i,{\theta }^{(n)})& =\frac{{\alpha }_{z_i}\varphi (x_i|{\theta }_{z_i}^{(n)})}{{\sum }_{k=1}^K{\alpha }_k\varphi (x_i|{\theta }_k^{(n)})}\\ Q(\theta ,{\theta }^{(n)})& =\sum _{i=1}^m[\sum _{z_i=j}^K\frac{{\alpha }_{z_i}\varphi (x_i|{\theta }_{z_i}^{(n)})}{{\sum }_{k=1}^K{\alpha }_k\varphi (x_i|{\theta }_k^{(n)})}\log {\alpha }_{z_i}\varphi (x_i|{\theta }_{z_i})]\end{array}$$

9.3.2书上的推导

下面的公式中，对于每个样本$x_j$的参数$r_{j1},r_{j2},\cdots ,r_{jK}$中，只有一个值为1，其他值都为0。
那么每个样本的概率可以表示为
$$p(x_j,r_j|\theta )=\prod _{k=1}^K[{\alpha }_k^{r_{jk}}\varphi (x_j|{\theta }_k{)}^{r_{jk}}]$$
所有样本的极大似然函数可以表示为
$$\begin{array}{rl}P(x,\gamma |\theta )& =\prod _{j=1}^N\prod _{k=1}^K[{\alpha }_k^{r_{jk}}\varphi (x_j|{\theta }_k{)}^{r_{jk}}]\\ \prod _{j=1}^N\prod _{k=1}^K{\alpha }_k^{r_{jk}}& =\prod _{k=1}^K\prod _{j=1}^N{\alpha }_k^{r_{jk}}=\prod _{k=1}^K{\alpha }_k^{{\sum }_{j=1}^N{r}_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\text{令nk为依赖第k个模型生成观测值的样本数}\\ P(x,\gamma |\theta )& =\prod _{k=1}^K[{\alpha }_k^{n_k}\prod _{j=1}^N\varphi (x_j|{\theta }_k{)}^{r_{jk}}]\\ & =\prod _{k=1}^K[{\alpha }_k^{n_k}\prod _{j=1}^N[\frac{1}{\sqrt{2\pi }{\sigma }_k}exp(-\frac{(x_j-{\mu }_k{)}^2}{2{\sigma }_k^2}){]}^{r_{jk}}]\end{array}$$
其中$n_k={\sum }_{j=1}^N{r}_{jk},{\sum }_{k=1}^K{n}_k=N$

对数似然函数
$$\begin{array}{r}\log P(x,\gamma |\theta )=\sum _{k=1}^K\{n_k\log {\alpha }_k+\sum _{j=1}^N{r}_{jk}[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(x_j-{\mu }_k{)}^2]\}\end{array}$$
算法E步
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E_{\gamma }[\log P(x,\gamma |\theta )|x,{\theta }^{(i)}]\\ & =E_{\gamma }\left\{\sum _{k=1}^K\right[n_k\log {\alpha }_k+\sum _{j=1}^N{r}_{jk}[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(x_j-{\mu }_k{)}^2]]\}\\ & =\sum _{k=1}^K[E(n_k)\log {\alpha }_k+\sum _{j=1}^{N}E(r_{jk})[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(x_j-{\mu }_k{)}^2]]\end{array}$$
计算
$$\begin{array}{rl}E(n_k)& =E(\sum _{j=1}^N{r}_{jk})=\sum _{j=1}^{N}E(r_{jk})\\ E(r_{jk})& =E(r_{jk}|x,{\theta }^{(i)})=P(r_{jk}=1|x,{\theta }^{(i)})=\frac{P(r_{jk}=1,x|{\theta }^{(i)})}{P(x|{\theta }^{(i)})}\\ & =\frac{P(r_{jk}=1,x|{\theta }^{(i)})}{{\sum }_{k=1}^{K}P(r_{jk}=1,x|{\theta }^{(i)})}\\ & =\frac{P(r_{jk}=1|{\theta }^{(i)})P(x|r_{jk}=1,{\theta }^{(i)})}{{\sum }_{k=1}^{K}P(r_{jk}=1|{\theta }^{(i)})P(x|r_{jk}=1,{\theta }^{(i)})}\\ & =\frac{{\alpha }_k^{(i)}\varphi (x_j,{\theta }_k^{(i)})}{{\sum }_{k=1}^K{\alpha }_k^{(i)}\varphi (x_j,{\theta }_k^{(i)})}\end{array}$$
带入Q函数得到
$$Q(\theta ,{\theta }^{(i)})=\sum _{k=1}^K[\sum _{j=1}^{N}E(r_{jk})\log {\alpha }_k+\sum _{j=1}^{N}E(r_{jk})[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(x_j-{\mu }_k{)}^2]]$$

---

两种方法结果等价
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =\sum _{k=1}^K[\sum _{j=1}^{N}E(r_{jk})\log {\alpha }_k+\sum _{j=1}^{N}E(r_{jk})[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(x_j-{\mu }_k{)}^2]]\\ & =\sum _{k=1}^K[\sum _{j=1}^{N}E(r_{jk})\log {\alpha }_k\varphi (x_j|\theta )]\\ & =\sum _{j=1}^N[\sum _{k=1}^K\frac{{\alpha }_k^{(i)}\varphi (x_j,{\theta }_k^{(i)})}{{\sum }_{l=1}^K{\alpha }_l^{(i)}\varphi (x_j,{\theta }_l^{(i)})}\log {\alpha }_k\varphi (x_j|\theta )]\end{array}$$

---

算法M步

• 对于${\mu }_k,{\sigma }_k^2$，没有约束条件，直接对Q函数求导
  $$\begin{array}{rl}\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\mu }_k}& =\sum _{j=1}^{N}E(r_{jk})[\frac{1}{2{\sigma }^2}(x_j-{\mu }_k)]=0\\ \frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\sigma }_k^2}& =\sum _{j=1}^{N}E(r_{jk})[-\frac{1}{2{\sigma }_k^2}+\frac{1}{2{\sigma }_k^4}(x_j-{\mu }_k{)}^2]=0\end{array}$$
  得到
  $$\begin{array}{rl}{\mu }_k& =\frac{{\sum }_{j=1}^{N}E(r_{jk})x_j}{{\sum }_{j=1}^{N}E(r_{jk})}\\ {\sigma }_k& =\frac{{\sum }_{j=1}^{N}E(r_{jk})(x_j-{\mu }_k{)}^2}{{\sum }_{j=1}^{N}E(r_{jk})}\end{array}$$
• 对于${\alpha }_k$,由于存在条件${\sum }_{k=1}^K{\alpha }_k=1$,原问题的拉格朗日函数为
  $$\begin{array}{rl}F(\alpha )& =Q(\alpha )+\beta (\sum _{k=1}^K{\alpha }_k-1)=\sum _{k=1}^K\sum _{j=1}^{N}E(r_{jk})\log {\alpha }_k+G()+\beta (\sum _{k=1}^K{\alpha }_k-1)\\ \frac{\partial F(\alpha )}{\partial {\alpha }_k}& =\frac{{\sum }_{j=1}^{N}E(r_{jk})}{{\alpha }_k}+\beta =0\to {\alpha }_k=\frac{{\sum }_{j=1}^{N}E(r_{jk})}{\beta }\\ \sum _{k=1}^K{\alpha }_k& =\frac{{\sum }_{k=1}^K{\sum }_{j=1}^{N}E(r_{jk})}{\beta }=1\to \beta =\sum _{k=1}^K\sum _{j=1}^{N}E(r_{jk})=N\end{array}$$
  得到
  $${\alpha }_k^2=\frac{{\sum }_{j=1}^{N}E(r_{jk})}{N}$$

9.3.3 西瓜书的推导

西瓜书上并没有利用Q函数，直接拉格朗日函数求导了
$$\begin{array}{rl}P(x_i|\theta )& =\sum _{k=1}^K{\alpha }_k\varphi (x_i|{\theta }_k)\\ L& =\sum _{j=1}^N\log \sum _{k=1}^K{\alpha }_k\varphi (x_i|{\theta }_k)\\ s.t.& \sum _{k=1}^K{\alpha }_k=1\end{array}$$
拉格朗日函数为
$$l=\sum _{j=1}^N\log \sum _{k=1}^K{\alpha }_k\varphi (x_i|{\theta }_k)+\lambda (\sum _{k=1}^K{\alpha }_k-1)$$
对${\alpha }_k,{\mu }_k,{\Sigma }_k$求导会和用EM算法得到一样的结果

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

输入：观测数据$x_i,\cdots ,x_N$，高斯混合模型
输出：高斯混合模型参数

• 1)取参数的初始值开始迭代
• 2)E步：依据当前模型参数，计算分模型k对观测数据的影响度
  $${\hat{\gamma }}_{jk}=E({\gamma }_{jk})=\frac{{\alpha }_k^{(i)}\varphi (x_j,{\theta }_k^{(i)})}{{\sum }_{k=1}^K{\alpha }_k^{(i)}\varphi (x_j,{\theta }_k^{(i)})},j=1,2,\cdots ,N;k=1,2,\cdots ,K$$
• 3)M步：计算新一轮的模型参数
  $${\mu }_k=\frac{{\sum }_{j=1}^{N}E(r_{jk})x_j}{{\sum }_{j=1}^{N}E(r_{jk})}$$
  $${\sigma }_k^2=\frac{{\sum }_{j=1}^{N}E(r_{jk})(x_j-{\mu }_k{)}^2}{{\sum }_{j=1}^{N}E(r_{jk})}$$
  $${\alpha }_k=\frac{{\sum }_{j=1}^{N}E(r_{jk})}{N}$$
• 4)重复2,3步至收敛