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

文章目录

9.EM算法及其推广9.1概念9.1.1模型解释9.1.2极大似然函数 9.2EM算法9.2.1推导9.2.2算法9.2.3 算法的收敛性 9.3EM在高斯混合模型中的应用9.3.1利用EM结论的推导9.3.2书上的推导9.3.3 西瓜书的推导9.3.4高斯混合模型参数估计的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步至收敛