机器学习（九）——EM算法
http://antkillerfarm.github.io/
高斯混合模型和EM算法（续）
E-Step中，根据贝叶斯公式可得：
p(z^{(i)}=j\vert x^{(i)};\phi,\mu,\Sigma)=\frac{p(x^{(i)}\vert z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{\sum_{l=1}^kp(x^{(i)}\vert z^{(i)}=l;\mu,\Sigma)p(z^{(i)}=l;\phi)}
将假设模型的各概率密度函数代入上式，即可计算得到w(i)jw_j^{(i)}。
相比于K-means算法，GMM算法中的z(i)z^{(i)}是个概率值，而非确定值，因此也被称为soft assignments。
K-means算法各个聚类的特征都是一样的，也就是“圆圈”的半径一致。而GMM算法的“圆圈”半径可以不同。如下面两图所示：
|
|
|---|---|
| K-means算法 | GMM算法 |
GMM算法结果也是局部最优解。对其他参数取不同的初始值进行多次计算同样适用于GMM算法。
参考：
http://cseweb.ucsd.edu/~atsmith/project1_253.pdf
EM算法
本节将进一步讨论EM算法的性质，并将之应用到使用latent random variables的一大类估计问题中。
Jensen不等式
首先我们给出凸函数的定义：
设f是定义域为实数的函数，如果对于所有的实数x，f′′(x)≥0f''(x)\ge 0，那么f是凸函数。当x是向量时，如果其Hessian矩阵H是半正定的（H≥0H\ge 0），那么f是凸函数。如果f′′(x)>0f''(x)>0或H>0H>0，那么称f是严格凸函数。
Jensen不等式表述如下:
如果f是凸函数，X是随机变量，那么E[f(X)]≥f(EX)E[f(X)]\ge f(EX)。
特别地，如果f是严格凸函数，那么E[f(X)]=f(EX)E[f(X)]=f(EX)，当且仅当p(X=EX)=1p(X=EX)=1。也就是说X是常量。
在不引起误会的情况下，E[X]E[X]简写做EXEX。
这个不等式的含义如下图所示：
Jensen不等式应用于凹函数时，不等号方向反向，也就是E[f(X)]≤f(EX)E[f(X)]\le f(EX)
注：Johan Ludwig William Valdemar Jensen，1859～1925，丹麦人。主业工程师，副业数学家。
EM算法的一般形式
EM算法一般形式的似然函数为：
\ell(\theta)=\sum_{i=1}^m\log p(x^{(i)};\theta)=\sum_{i=1}^m\log \sum_zp(x^{(i)},z^{(i)};\theta)\tag{1}
根据这个公式直接求θ\theta一般比较困难，但是如果确定了z之后，求解就容易了。
EM算法是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化ℓ(θ)\ell(\theta)，我们可以不断地建立ℓ(θ)\ell(\theta)的下界(E-Step),然后优化下界(M-Step)。
这里首先假设z(i)z^{(i)}的分布为QiQ_i。显然：
\sum_zQ_i(z)=1,Q_i(z)\ge 0
\begin{align} \ell(\theta)&=\sum_{i=1}^m\log \sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta) \\&=\sum_{i=1}^m\log \sum_{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \\&\ge\sum_{i=1}^m\sum_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\tag{2} \end{align}
这里解释一下最后一步的不等式变换。
首先，根据数学期望的定义公式：
E[X]=\sum_{i=1}^nx_ip_i
可知：
\sum_{z^{(i)}}Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=E\left[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\right]
又因为f(x)=logxf(x)=\log x是凹函数，根据Jensen不等式可得：
f\left(E_{z^{(i)}\sim Q_i}\left[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\right]\right)\ge E_{z^{(i)}\sim Q_i}\left[f\left(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\right)\right]
综上，公式2给出了ℓ(θ)\ell(\theta)的下界。对于QiQ_i的选择，有多种可能，那种更好的？
假设θ\theta已经给定，那么ℓ(θ)\ell(\theta)的值就取决于Qi(z(i))Q_i(z^{(i)})和p(x(i),z(i))p(x^{(i)},z^{(i)})了。我们可以通过调整这两个概率，使下界不断上升，以逼近ℓ(θ)\ell(\theta)的真实值。当不等式变成等式时，下界达到最大值。
由Jensen不等式相等的条件可知，ℓ(θ)\ell(\theta)的下界达到最大值的条件为：
\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c
其中c为常数。
这实际上表明：
Q_i(z^{(i)})\propto p(x^{(i)},z^{(i)};\theta)
其中的∝\propto符号是两者成正比例的意思。
从中还可以推导出：
\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c=\frac{\sum_zp(x^{(i)},z^{(i)};\theta)}{\sum_zQ_i(z^{(i)})}
因为∑zQi(z(i))=1\sum_zQ_i(z^{(i)})=1，所以上式可变形为：
Q_i(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta)}{\sum_zp(x^{(i)},z^{(i)};\theta)}=\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)}=p(z^{(i)}\vert x^{(i)};\theta)
可见，当Qi(z(i))Q_i(z^{(i)})为z(i)z^{(i)}的后验分布时，ℓ(θ)\ell(\theta)的下界达到最大值。
因此，EM算法的过程为：
Repeat until convergence {
(E-step) For each i：
Qi(z(i)):=p(z(i)|x(i);θ)Q_i(z^{(i)}):=p(z^{(i)}\vert x^{(i)};\theta)
(M-step) Update the parameters：
θ:=argmaxθ∑i∑zQi(z(i))logp(x(i),z(i);θ)Qi(z(i))\theta:=\arg\max_\theta\sum_i\sum_zQ_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}
}
如何保证算法的收敛性呢？
如果我们用θ(t)\theta^{(t)}和θ(t+1)\theta^{(t+1)}表示EM算法第t次和第t+1次迭代后的结果，那么我们的任务就是证明ℓ(θ(t))≤ℓ(θ(t+1))\ell(\theta^{(t)})\le \ell(\theta^{(t+1)})。
由公式1和2可得：
\ell(\theta)\ge\sum_i\sum_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\tag{3}
由之前的讨论可以看出，E-Step中的步骤是使上式的等号成立的条件，即：
\ell(\theta^{(t)})=\sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})}
因为公式3对于任意QiQ_i和θ\theta都成立，因此：
\ell(\theta^{(t+1)})\ge\sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})}
因为M-Step的最大化过程，可得：
\sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})}\ge\sum_i\sum_{z^{(i)}}Q_i^{(t)}(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})}
综上可得：ℓ(θ(t))≤ℓ(θ(t+1))\ell(\theta^{(t)})\le \ell(\theta^{(t+1)})
事实上，如果我们定义：
J(Q,\theta)=\sum_i\sum_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}
则EM算法可以看作是J函数的坐标上升法。E-Step固定θ\theta，优化Q；M-Step固定Q，优化θ\theta。
重新审视混合高斯模型
下面给出混合高斯模型各参数的推导过程。
E-Step很简单：
w_j^{(i)}=Q_i(z^{(i)}=j)=p(z^{(i)}=j\vert x^{(i)};\phi,\mu,\Sigma)
在M-Step中，我们将各个变量和分布的概率密度函数代入，可得：
\begin{align} &\sum_{i=1}^m\sum_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \\&=\sum_{i=1}^m\sum_{j=1}^kQ_i(z^{(i)}=j)\log\frac{p(x^{(i)}\vert z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{Q_i(z^{(i)}=j)} \\&=\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}\log\frac{\frac{1}{(2\pi)^{n/2}\lvert\Sigma_j\rvert^{1/2}}\exp\left(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)\right)\cdot \phi_j}{w_j^{(i)}} \end{align}
对μl\mu_l求导，可得：
\begin{align} &\nabla_{\mu_l}\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}\log\frac{\frac{1}{(2\pi)^{n/2}\lvert\Sigma_j\rvert^{1/2}}\exp\left(-\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j)\right)\cdot \phi_j}{w_j^{(i)}} \\&=-\nabla_{\mu_l}\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}\frac{1}{2}(x^{(i)}-\mu_j)^T\Sigma_j^{-1}(x^{(i)}-\mu_j) \\&=\frac{1}{2}\sum_{i=1}^mw_l^{(i)}\nabla_{\mu_l}(2\mu_l^T\Sigma_l^{-1}x^{(i)}-\mu_l^T\Sigma_l^{-1}\mu_l)=\sum_{i=1}^mw_l^{(i)}(\Sigma_l^{-1}x^{(i)}-\Sigma_l^{-1}\mu_l)\tag{4} \end{align}
这里对最后一步的推导，做一个说明。为了便于以下的讨论，我们引入符号“tr”，该符号表示矩阵的主对角线元素之和，也被叫做矩阵的“迹”（Trace）。按照通常的写法，在不至于误会的情况下，tr后面的括号会被省略。
tr的其他性质还包括：
\operatorname{tr}\,a=a,a\in R\tag{5.1}
\operatorname{tr}(A + B) = \operatorname{tr}(A) + \operatorname{tr}(B)\tag{5.2}
\operatorname{tr}(cA) = c \operatorname{tr}(A)\tag{5.3}
\operatorname{tr}(A) = \operatorname{tr}(A^{\mathrm T})\tag{5.4}
\operatorname{tr}(AB) = \operatorname{tr}(BA)\tag{5.5}
\operatorname{tr}(ABCD) = \operatorname{tr}(BCDA) = \operatorname{tr}(CDAB) = \operatorname{tr}(DABC)\tag{5.6}
\nabla_A\operatorname{tr}(AB)=B^T\tag{5.7}
\nabla_{A^T}f(A)=(\nabla_Af(A))^T\tag{5.8}
\nabla_A\operatorname{tr}(ABA^TC)=CAB+C^TAB^T\tag{5.9}
\nabla_{A^T}\operatorname{tr}(ABA^TC)=B^TA^TC^T+BA^TC\tag{5.10}
\nabla_A|A|=|A|(A^{-1})^T\tag{5.11}
因为μTlΣ−1lμl\mu_l^T\Sigma_l^{-1}\mu_l是实数，由公式5.1可得：
\nabla_{\mu_l}\mu_l^T\Sigma_l^{-1}\mu_l=\nabla_{\mu_l}\operatorname{tr}(\mu_l^T\Sigma_l^{-1}\mu_l)
由公式5.10可得：
\nabla_{\mu_l}\operatorname{tr}(\mu_l^T\Sigma_l^{-1}\mu_l)=(\Sigma_l^{-1})^T(\mu_l^T)^TI^T+\Sigma_l^{-1}(\mu_l^T)^TI=(\Sigma_l^{-1})^T\mu_l+\Sigma_l^{-1}\mu_l
因为Σ−1l\Sigma_l^{-1}是对称矩阵，因此，综上可得：
\nabla_{\mu_l}\mu_l^T\Sigma_l^{-1}\mu_l=2\Sigma_l^{-1}\mu_l\tag{5.12}
回到正题，令公式4等于0，可得：
\mu_j:=\frac{\sum_{i=1}^mw_j^{(i)}x^{(i)}}{\sum_{i=1}^mw_j^{(i)}}
同样的，对ϕj\phi_j求导，可得：
\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}\log\phi_j
因为存在∑kj=1ϕj=1\sum_{j=1}^k\phi_j=1这样的约束，因此需要使用拉格朗日函数：
\mathcal{L}(\phi)=\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}\log\phi_j+\beta(\sum_{j=1}^k\phi_j-1)
这里的β\beta是拉格朗日乘子，ϕj≥0\phi_j\ge 0的约束条件不用考虑，因为对logϕj\log\phi_j求导已经隐含了这个条件。
因此：
\frac{\partial\mathcal{L}(\phi)}{\partial\phi_j}=\sum_{i=1}^m\frac{w_j^{(i)}}{\phi_j}+\beta
令上式等于0，可得：
\frac{\sum_{i=1}^mw_j^{(i)}}{\phi_j}=-\beta=\frac{\sum_{i=1}^m\sum_{j=1}^kw_j^{(i)}}{\sum_{j=1}^k\phi_j}
因为∑kj=1ϕj=1,∑kj=1w(i)j=1\sum_{j=1}^k\phi_j=1,\sum_{j=1}^kw_j^{(i)}=1，所以：
-\beta=\frac{\sum_{i=1}^m1}{1}=m
因此：
\phi_j:=\frac{1}{m}\sum_{i=1}^mw_j^{(i)}
因子分析
之前的讨论都是基于样本个数m远大于特征数n的，现在来看看m≪nm\ll n的情况。
这种情况本质上意味着，样本只覆盖了很小一部分的特征空间。当我们应用高斯模型的时候，会发现协方差矩阵Σ\Sigma根本就不存在，自然也就没法利用之前的方法了。
那么我们应该怎么办呢？
对Σ\Sigma的限制
有两种方法可以对Σ\Sigma进行限制。
方法一：
设定Σ\Sigma为对角线矩阵，即所有非对角线元素都是0。其对角线元素为：
\Sigma_{jj}=\frac{1}{m}\sum_{i=1}^m(x_j^{(i)}-\mu_j)^2
二维多元高斯分布在平面上的投影是个椭圆，中心点由μ\mu决定，椭圆的形状由Σ\Sigma决定。Σ\Sigma如果变成对角阵，就意味着椭圆的两个轴都和坐标轴平行了。
方法二：
还可以进一步约束Σ\Sigma，可以假设对角线上的元素都是相等的，即：
\Sigma=\sigma^2I
其中：
\sigma=\frac{1}{mn}\sum_{j=1}^n\sum_{i=1}^m(x_j^{(i)}-\mu_j)^2
机器学习（九）——EM算法相关推荐
- 机器学习九大算法---支持向量机
机器学习九大算法---支持向量机 出处:结构之法算法之道blog. 前言 动笔写这个支持向量机(support vector machine)是费了不少劲和困难的,原因很简单,一者这个东西本身就并不好 ...
- em算法 实例 正态分布_【机器学习】EM算法详细推导和讲解
今天不太想学习,炒个冷饭,讲讲机器学习十大算法里有名的EM算法,文章里面有些个人理解,如有错漏,还请读者不吝赐教. 众所周知,极大似然估计是一种应用很广泛的参数估计方法.例如我手头有一些东北人的身高的 ...
- em算法怎么对应原有分类_机器学习基础-EM算法
EM算法也称期望最大化(Expectation-Maximum,简称EM)算法,它是一个基础算法,是很多机器学习领域算法的基础,比如隐式马尔科夫算法(HMM), LDA主题模型的变分推断等等.本文就对 ...
- 机器学习算法_机器学习之EM算法和概率图模型
[晓白]今天我准备更新Machine Learning系列文章希望对机器学习复习和准备面试的同学有帮助!之前更新了感知机和SVM,决策树&代码实战,关注我的专栏可以的文章哦!今天继续更新EM算 ...
- 机器学习基础 EM算法
文章目录 一.初识EM算法 二.EM算法介绍 1. 极大似然估计 1.1 问题描述 1.2 用数学知识解决现实问题 1.3 最大似然函数估计值的求解步骤 2. EM算法实例描述 3. EM算法流程 三 ...
- 机器学习 之 EM算法
EM(Expectation-Maximum)算法也称期望最大化算法 EM算法是最常见的隐变量估计方法,在机器学习中有极为广泛的用途,例如常被用来学习高斯混合模型(Gaussian mixture m ...
- 机器学习之EM算法的原理推导及相关知识总结
文章目录 1.知道先验概率和后验概率 2.了解高斯混合模型GMM 3.通过最大似然估计推导EM算法的过程的实例 4.EM算法 5.知道pLSA模型 1.知道先验概率和后验概率 先验概率(prior p ...
- 【机器学习】EM算法
EM算法 目录 一.似然函数与极大似然估计 二.Jenson不等式 三.数学期望的相关定理 四.边缘分布列 五.坐标上升法 六.EM算法 1. 概论 2. 算法流程 3. 算法的推导 4. 敛散性证明 ...
- 机器学习之EM算法的原理及推导(三硬币模型)及Python实现
EM算法的简介 EM算法由两步组成:E步和M步,是最常用的迭代算法. 本文主要参考了李航博士的<统计学习方法> 在此基础上主要依据EM算法原理补充了三硬币模型的推导. 1.EM算法的原理 ...
- 【机器学习】EM 算法
在本文中规定,下凸函数为类似于 f ( x ) = x 2 f(x)=x^2 f(x)=
最新文章
- 因女朋友的一个建议，这位程序员创立仅 551 天公司就被 10 亿美元收购了
- 什么？欧洲也有个恩智浦杯?
- 7. Docker - 网络管理
- 7 项目人力资源管理
- 详解如何充分发挥先验信息优势，用MRC框架解决各类NLP任务
- 控制SAP Spartacus shipping address页面spinner显示的逻辑
- cocos2d-x学习笔记之图片分辨率适配
- LeetCode算法题-Delete Node in a Linked List（Java实现）
- 四选1数据选择器的VHDL程序设计
- PB系统连接原生数据库DB文件
- 计算机word表格三线形,word中制作三线表格的四种方法
- 仿微信发送位置(高德地图定位)
- AIT Worldwide Logistics与卡利塔航空建立民航预备队合作伙伴关系
- 关于App的数据统计和分析
- ctf你好 1.1:灵魂拷问
- react native 啧啧啧
- 如何安装numpy库
- Excel2016设置下拉选项并自动匹配单元格颜色
- 小度路由器离线下载根本就是垃圾
- VASP输入INCAR文件