文章目录

EM 算法
1. EM 算法的引入
- 三硬币模型
2. EM 算法
- Q 函数
参考文献

机器学习 | 目录

无监督学习 | GMM 高斯混合原理及Sklearn实现

本文大部分内容搬运自李航老师的《统计学习方法》^[1]，以给出 EM 算法较为完整的定义。

EM 算法

EM 算法是一种迭代算法，1977 年由 Dempster 等人总结提出，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计，或极大后验估计。

EM 算法的每次迭代由两步组成：E 步：求期望（expectation）；M 步：求极大（maximization）。所以这一算法称为期望极大算法（expectation maximization algorithm, EM）。

1. EM 算法的引入

概率模型有时即含有观测数据（observable varible，已知），又含有隐变量或潜在变量（latent varible，未知），如果概率模型的变量都是观测变量，那么给定数据，可以直接使用极大似然估计，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单地使用这些估计方法。EM 算法就是含有隐变量的概率模型参数的极大似然估计。

三硬币模型

假设有三枚硬币，分别记作 A，B，C。这些硬币正面出现的概率分别是 π\piπ，ppp 和 qqq。进行如下掷硬币试验：

先掷硬币 A ，若为正面则选硬币 B ，若为反面则选硬币 C ；然后掷 A 选出的硬币（ B 或 C ），若为正面则记作 1，若为反面则记作 0，试验结束。

独立重复试验 nnn 次（这里 n=10n=10n=10），观测结果如下：

1,1,0,1,0,0,1,0,1,11,1,0,1,0,0,1,0,1,11,1,0,1,0,0,1,0,1,1

假设只能观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三硬币正面出现的概率，即三硬币模型的参数。

解三硬币模型可以写作：

P(y∣θ)=∑zP(y,z∣θ)=∑zP(z∣θ)P(y∣z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y(1)\begin{aligned} P(y | \theta) &=\sum_{z} P(y, z | \theta)=\sum_{z} P(z | \theta) P(y | z, \theta) \\ &=\pi p^{y}(1-p)^{1-y}+(1-\pi) q^{y}(1-q)^{1-y} \end{aligned} \tag{1}P(y∣θ)=z∑P(y,z∣θ)=z∑P(z∣θ)P(y∣z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y(1)

这里，随机变量 yyy 是观测变量，表示一次试验观测的结果是 1 或 0；随机变量 zzz 是隐变量，表示未观测到的掷硬币 A 的结果；θ=(π,p,q)\theta=(\pi,p,q)θ=(π,p,q) 是模型参数。这一模型是以上数据的生成模型。注意，随机变量 yyy 的数据可以观测，随机变量 zzz 的数据不可观测。

将观测数据表示为 Y=(Y1,Y2,...,Yn)TY=(Y_1,Y_2,...,Y_n)^TY=(Y1,Y2,...,Yn)T ，未观测数据表示为 Z=(Z1,Z2,...,Zn)TZ=(Z_1,Z_2,...,Z_n)^TZ=(Z1,Z2,...,Zn)T，则观测数据的似然函数为：

P(Y∣θ)=∑ZP(Z∣θ)P(Y∣Z,θ)(2)P(Y|\theta)=\sum_Z P(Z|\theta)P(Y|Z,\theta) \tag{2}P(Y∣θ)=Z∑P(Z∣θ)P(Y∣Z,θ)(2)

即

P(Y∣θ)=∏j=1n[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj](3)P(Y | \theta)=\prod_{j=1}^{n}\left[\pi p^{y_{j}}(1-p)^{1-y_{j}}+(1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}\right] \tag{3} P(Y∣θ)=j=1∏n[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj](3)

考虑求模型参数 θ=(π,p,q)\theta=(\pi,p,q)θ=(π,p,q)，即：

θ^=argmax⁡θlogP(Y∣θ)(4)\hat{\theta}=arg \max \limits_{\theta} logP(Y|\theta) \tag{4}θ^=argθmaxlogP(Y∣θ)(4)

这个问题没有解析解，只有通过迭代的方法求解。EM 算法就是可以用于求解这个问题的一种迭代算法。下面给出针对以上问题的 EM 算法，其推导过程省略。

EM 算法首先选取参数的初值，记作 θ(0)=(π(0),p(0),q(0))\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})θ(0)=(π(0),p(0),q(0)) ，然后通过下面的步骤迭代计算参数的估计值，直至收敛为止。第 iii 次迭代参数的估计值为 θ(i)=(π(i),p(i),q(i))\theta^{(i)}=(\pi^{(i)},p^{(i)},q^{(i)})θ(i)=(π(i),p(i),q(i))。EM 算法的第 i+1i+1i+1 次迭代如下：

E 步：计算在模型参数 π(i),p(i),q(i)\pi^{(i)},p^{(i)},q^{(i)}π(i),p(i),q(i) 下观测数据 yiy_iyi 来自掷硬币 B 的概率：

μ(i+1)=π(i)(p(i))yj(1−p(i))1−yjπ(i)(p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))1−yj(5)\mu^{(i+1)}=\frac{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}}{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}+\left(1-\pi^{(i)}\right)\left(q^{(i)}\right)^{y_{j}}\left(1-q^{(i)}\right)^{1-y_{j}}} \tag{5}μ(i+1)=π(i)(p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))1−yjπ(i)(p(i))yj(1−p(i))1−yj(5)

M 步：计算模型参数的新估计值

π(i+1)=1n∑j=1nμj(i+1)(6)\pi^{(i+1)}=\frac{1}{n} \sum_{j=1}^{n} \mu_{j}^{(i+1)} \tag{6}π(i+1)=n1j=1∑nμj(i+1)(6)

p(i+1)=∑j=1nμj(i+1)yj∑j=1nμj(i+1)(7)p^{(i+1)}=\frac{\sum_{j=1}^{n} \mu_{j}^{(i+1)} y_{j}}{\sum_{j=1}^{n} \mu_{j}^{(i+1)}} \tag{7}p(i+1)=∑j=1nμj(i+1)∑j=1nμj(i+1)yj(7)

q(i+1)=∑j=1n(1−μj(i+1))yj∑j=1n(8)q^{(i+1)}=\frac{\sum_{j=1}^{n}\left(1-\mu_{j}^{(i+1)}\right) y_{j}}{\sum_{j=1}^{n}} \tag{8}q(i+1)=∑j=1n∑j=1n(1−μj(i+1))yj(8)

进行数字计算。假设模型参数的初值取为：

π(0)=0.5,p(0)=0.5,q(0)=0.5\pi^{(0)}=0.5,\quad p^{(0)}=0.5,\quad q^{(0)}=0.5π(0)=0.5,p(0)=0.5,q(0)=0.5

由式 (5)，对 yj=1y_j=1yj=1 与 yj=0y_j=0yj=0 均有 μj(1)=0.5\mu_j^{(1)}=0.5μj(1)=0.5。

利用迭代公式 (6-8) ，得到：

π(1)=0.5,p(1)=0.6,q(1)=0.6\pi^{(1)}=0.5,\quad p^{(1)}=0.6,\quad q^{(1)}=0.6π(1)=0.5,p(1)=0.6,q(1)=0.6

由式 (5)，

μj(2)=0.5,j=1,2,...,10\mu_j^{(2)}=0.5, \quad j=1,2,...,10μj(2)=0.5,j=1,2,...,10

继续迭代，得：

π(2)=0.5,p(2)=0.6,q(2)=0.6\pi^{(2)}=0.5,\quad p^{(2)}=0.6,\quad q^{(2)}=0.6π(2)=0.5,p(2)=0.6,q(2)=0.6

可以看到参数以及收敛，于是得到模型参数 θ\thetaθ 的极大似然估计：

π^=0.5,p^=0.6q^=0.6\hat{\pi}=0.5, \quad \hat{p}=0.6 \quad \hat{q}=0.6π^=0.5,p^=0.6q^=0.6

π=0.5\pi = 0.5π=0.5 表示硬币 A 是均匀的，这一结果容易理解。

如果取初值π(0)=0.5,p(0)=0.6,q(0)=0.7\pi^{(0)}=0.5,\quad p^{(0)}=0.6,\quad q^{(0)}=0.7π(0)=0.5,p(0)=0.6,q(0)=0.7，那么得到的模型参数的极大似然估计是 π^=0.4064,p^=0.5368q^=0.6432\hat{\pi}=0.4064, \quad \hat{p}=0.5368 \quad \hat{q}=0.6432π^=0.4064,p^=0.5368q^=0.6432 。这就是说，EM 算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

一般地，用 YYY 表示观测随机变量的数据，ZZZ 表示隐随机变量的数据。YYY 和 ZZZ 连在一起称为完全数据（complete-data），观测数据 YYY 又称为不完全数据（incomplete-data）。假设给定观测数据 YYY，其概率分布是 P(Y∣θ)P(Y|\theta)P(Y∣θ)，其中 θ\thetaθ 是需要估计的模型参数，那么不完全数据 YYY 的似然函数是 P(Y∣θ)P(Y|\theta)P(Y∣θ)，对数似然函数为 L(θ)=log⁡P(Y∣θ)L(\theta)=\log P(Y|\theta)L(θ)=logP(Y∣θ)；假设 YYY 和 ZZZ 的联合概率分布是 P(Y,Z∣θ)P(Y,Z|\theta)P(Y,Z∣θ)，那么完全数据的对数似然函数是 log⁡P(Y,Z∣θ)\log P(Y,Z|\theta)logP(Y,Z∣θ)。

EM 算法通过迭代求 L(θ)=log⁡P(Y∣θ)L(\theta)=\log P(Y|\theta)L(θ)=logP(Y∣θ) 的极大似然估计。每次迭代包含两步：E 步，求期望；M 步，求极大化。

2. EM 算法

输入：观测变量数据 YYY，隐变量数据 ZZZ，联合分布 P(Y,Z∣θ)P(Y,Z|\theta)P(Y,Z∣θ)，条件分布 P(Z∣Y.θ)P(Z|Y.\theta)P(Z∣Y.θ)；

输出；模型参数 θ\thetaθ。

（1）选择参数的初值 θ(0)\theta^{(0)}θ(0)，开始迭代；

（2）E 步：记 θ(i)\theta^{(i)}θ(i) 为第 iii 次迭代参数 θ\thetaθ 的估计值，在第 i+1i+1i+1 次迭代的 E 步，计算：

Q(θ,θ(i))=EZ[log⁡P(Y,Z∣θ)∣Y,θ(i)]=∑Zlog⁡P(Y,Z∣θ)P(Z∣Y,θ(i))(9)\begin{aligned} Q\left(\theta, \theta^{(i)}\right) &=E_{Z}\left[\log P(Y, Z | \theta) | Y, \theta^{(i)}\right] \\ &=\sum_{Z} \log P(Y, Z | \theta) P\left(Z | Y, \theta^{(i)}\right) \end{aligned} \tag{9} Q(θ,θ(i))=EZ[logP(Y,Z∣θ)∣Y,θ(i)]=Z∑logP(Y,Z∣θ)P(Z∣Y,θ(i))(9)

（3）M 步：求使 Q(θ,θ(i))Q(\theta,\theta^{(i)})Q(θ,θ(i)) 极大化的 θ\thetaθ，确定第 i+1i+1i+1 次迭代的参数的估计值 θ(i+1)\theta^{(i+1)}θ(i+1)：

θ(i+1)=argmax⁡θQ(θ,θ(i))(10)\theta^{(i+1)}=arg \max \limits_{\theta} Q(\theta,\theta^{(i)}) \tag{10}θ(i+1)=argθmaxQ(θ,θ(i))(10)

（4）重复第 (2-3) 步，直到收敛。

下面关于 EM 算法作几点说明：

步骤 (1) 参数的初值可以任意选择，但需注意 EM 算法对初值是敏感的；

步骤 (2) E 步求 Q(θ,θ(i))Q(\theta,\theta^{(i)})Q(θ,θ(i)) 。Q 函数式中 Z 是未观测数据，Y 是观测数据。注意，Q(θ,θ(i))Q(\theta,\theta^{(i)})Q(θ,θ(i)) 的第 1 个表示要极大化的参数，第 2 个变元表示参数的当前估计值。每次迭代实际在求 Q 函数及其极大；

步骤 (3) M 步求 Q(θ,θ(i))Q(\theta,\theta^{(i)})Q(θ,θ(i)) 的极大化，得到 θ(i+1)\theta^{(i+1)}θ(i+1)，完成一次迭代 θ(i)→θ(i+1)\theta^{(i)}\to \theta^{(i+1)}θ(i)→θ(i+1)，且每次迭代使似然函数增大或达到局部极值；

步骤 (4) 给出停止迭代的条件，一般是对较小的正值 ε1,ε2\varepsilon_1,\varepsilon_2ε1,ε2，若满足：

∥θ(i+1)−θ(i)∥<ε1or∥Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∥<ε2(11)\left\|\theta^{(i+1)}-\theta^{(i)}\right\|<\varepsilon_{1} \quad or \quad\left\|Q\left(\theta^{(i+1)}, \theta^{(i)}\right)-Q\left(\theta^{(i)}, \theta^{(i)}\right)\right\|<\varepsilon_{2}\tag{11} ∥∥∥θ(i+1)−θ(i)∥∥∥<ε1or∥∥∥Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∥∥∥<ε2(11)

则停止迭代。

式 (9) 的函数 Q(θ,θ(i))Q(\theta,\theta^{(i)})Q(θ,θ(i)) 是 EM 算法的核心，称为 Q 函数（Q function）。

Q 函数

完全数据的对数似然函数 log⁡P(Y,Z∣θ)\log P(Y,Z|\theta)logP(Y,Z∣θ) 关于在给定观测数据 YYY 和当前参数 θ(i)\theta^{(i)}θ(i) 下对未观测数据 Z 的条件概率分布 P(Z∣Y,θ(i))P(Z|Y,\theta^{(i)})P(Z∣Y,θ(i)) 的期望称为 Q 函数，即：

Q(θ,θ(i))=Ez[log⁡P(Y,Z∣θ)∣Y,θ(i)](11)Q\left(\theta, \theta^{(i)}\right)=E_{z}\left[\log P(Y, Z | \theta) | Y, \theta^{(i)}\right] \tag{11}Q(θ,θ(i))=Ez[logP(Y,Z∣θ)∣Y,θ(i)](11)

参考文献

[1] 周志华. 机器学习[M]. 北京: 清华大学出版社, 2016: 155-158.

机器学习 | EM 算法原理相关推荐

【机器学习】算法原理详细推导与实现(七):决策树算法
[机器学习]算法原理详细推导与实现(七):决策树算法在之前的文章中,对于介绍的分类算法有逻辑回归算法和朴素贝叶斯算法,这类算法都是二分类的分类器,但是往往只实际问题中yyy不仅仅只有{0,1}\{0 ...
机器学习 EM算法从极大似然法谈起
极大似然法 VS EM算法 EM(Expectation Maximization)算法是复杂化的极大似然法(Maximum Likelihood). 极大似然法问题中,所有观测值属于同一分布,需要估 ...
《机器学习：算法原理与编程实践》的读书笔记：SMO部分最难，大部分代码基于Scikit-Learn，决策树其实用处不大
机器学习:算法原理与编程实践目录 [隐藏] 1 机器学习的基础 2 中文文本分类 3 决策树的发展 4 推荐系统原理 5 梯度寻优 6 神经网络初步 7 预测的技术与哲学 8 万能分类器:SVM 9 ...
(转载)机器学习知识点(十四)EM算法原理
1.引言以前我们讨论的概率模型都是只含观测变量(observable variable), 即这些变量都是可以观测出来的,那么给定数据,可以直接使用极大似然估计的方法或者贝叶斯估计的方法:但是当模型 ...
高斯-赛得尔迭代式 c++_高斯混合模型(Gaussian Mixture Model)与EM算法原理(一)
高斯混合模型(Gaussian Mixture Model)是机器学习中一种常用的聚类算法,本文介绍了其原理,并推导了其参数估计的过程.主要参考Christopher M. Bishop的<Pa ...
机器学习 EM算法理解
文章目录前言极大似然问题描述参数估计极大似然估计极大似然应用极小化代价函数 EM算法 EM算法推导基础知识凸函数 Jensen不等式期望推导过程算法流程另一种呢理解 EM算法 ...
scikit-learn机器学习常用算法原理及编程实战（五）
决策树决策树是最经典的机器学习模型之一.预测速度快,可以处理类别型数据和连续型数据.通过本章读者可以掌握以下内容: 信息熵及信息增益的概念,以及决策树的分裂的原则: 决策树的创建及剪枝算法: sci ...
EM算法原理简析——图解
一. 扯淡转眼间毕业快一年了,这期间混了两份工作,从游戏开发到算法.感觉自己还是喜欢算法,可能是大学混了几年算法吧!所以不想浪费基础... 我是个懒得写博客的人,混了几年coding,写的博客不超过 ...
深入理解机器学习——EM算法/最大期望算法（Expectation-Maximization Algorithm, EM）
分类目录:<深入理解机器学习>总目录在前面的讨论中,我们一直假设训练样本所有属性变量的值都已被观测到,即训练样本是"完整"的.但在现实应用中往往会遇到"不完 ...

机器学习 | EM 算法原理