EM算法及其推广学习笔记
前言
在学习隐马尔科夫模型时,在学习算法中指出了Baum-Welch算法,来实现对隐马尔科夫模型参数的求解。在该学习算法中用到了EM算法,因此我们先来看看EM算法到底是何方神圣。可自己在学习EM算法时,又遇到了一个坑,什么是极大似然函数?因此,本文先介绍极大似然函数的相关概念,然后再对EM算法进行物理映射和实际数学推导。本文需要大量概率论知识,在数学推导关键处会贴出相关参考资料和博文,以备后期查询使用。
思考
- 什么是极大似然函数,用来解决什么实际问题?
- EM算法是什么,用来解决什么问题?
正文
极大似然估计
1.概念
在已知试验结果(即是样本)的情况下,用来估计满足这些样本分布的参数,把可能性最大的那个参数
θ
作为真实
θ∗
的参数估计。
2.实例
在学习极大似然函数时,一直在思考用什么样的例子才能把概念讲清楚。其实最大的困惑便在于似然估计是假设你有了一堆样本,你需要根据这些样本来猜出某些概率参数,从而使得这些样本在你面前的概率最大。用公式表示为:
这里显然做了一个最基本假设,即出现在我们面前的样本表征了最真实的物理世界,没有之一,且不考虑外界噪声对样本的干扰。现在,咱们举一个实际的例子,来理解似然估计的物理含义。
假设我们已知得肺癌的概率为总体的0.2,而不得肺癌的概率为0.8。在得肺癌中的人群中,我们又做了相应的统计,即抽烟人群占0.7,而不抽烟人群只有0.3。在没有肺癌的人群中,抽烟者占0.2,不抽烟者占0.8。然而很不幸,在真实世界中,我们遇到了一个倒霉蛋,他得了肺癌。问题来了,他是不是抽烟呢?
这个例子不算是真正的似然估计最贴切的实例,但基本思想可以拿来表征一下,方便下面对公式的理解。很显然,得肺癌背后的概率模型显然和抽烟相关性很大,即想要得肺癌的概率最大,模型必须拟合到 θ=0.7 ,使得获肺癌的概率最大。因此,我们可以大胆猜测,该肺癌患者是吸烟人群。详细的关于似然函数的物理含义可以参看博文从最大似然到EM算法浅解。
我们先给出数学定义,极大似然估计可以分为离散型和连续型两类。
3.离散型
设
x
为离散型随机变量,
4.连续型
设
x
为连续型随机变量,概率密度函数为
5.应用举例
继续来看例子,假设进行一个实验,实验次数为10次,每次实验成功率为0.2,那么不成功的概率为0.8,用
y
来表示成功的次数。由于前后的实验是相互独立的,所以可以计算得到成功的次数的概率密度为:
该式子分为两项因子,10次实验中有
y
次成功,那么即在10次中随意挑选
显然, f(y;θ) 是关于随机变量 y 和概率参数
当
θ=0.2
时,我们可以得到y取不同值的概率分布情况,如下图所示:
当
θ=0.7
时,我们可以得到y取不同值的概率分布情况,如下图所示:
好了,现在假设我们在实验室,开始完成某个实验,我们并不知道该实验成功的概率是多少,但做了10次实验后,我们只成功了2次,用高中的拿点概率知识拿来求解,那不就是实验成功率为0.2。的确,但由于实验次数相当的小,这里的0.2并非是真正的概率,而只是我们实验成功的频率。如抛一枚硬币,抛个10次,可能正面朝上的频率为0.6,但我们都知道,实际正面朝上的概率为0.5。那如何让频率接近0.5呢,不断的增加实验次数即可,你抛个2万次试试。所以我们不能简单的就把这个问题中求解的0.2作为我们的答案,我们也不可能大量重复实验来统计该实验成功率。遇到这种情况,我们便用到了似然估计方法。
如上给出了似然函数:
现在我们已知实验次数为2,我们要求 θ 使得该似然函数取到极大值,求极值问题只需要对 θ 求导即可,如果是多参量,那么可以对它求偏导,得到似然方程组,同样能求出想要的解。
求得 θ=0.2 。咦,算出来的答案是一样的,这不是多此一举嘛,但上述实验成功次数背后的参数 θ 模型是一维的,即我们可以用“肉眼”来直接看出答案,假如我们这次不是观察实验成功次数,来猜实验成功率,我换个问题问,假设我们班男生身高符合高斯分布,即男生身高概率密度函数符合高斯分布,给你一群男生的身高,请告诉我高斯分布的方差和均值分别是多少?这种情况下,其背后的 θ=[μ,π] 含有两个参数,简单的靠肉眼观察显然无法给出答案,因此我们需要借助数学工具,来理论化的证明说,当看到这一群男生的身高时,我们能找到参数 θ=[μ,π] 使得出现 这群男生身高的概率最大,注意是这群而不是那群!言外之意就是说, θ 参数能够对男生进行分类!隐约看出了EM算法中的一些思想。
我们再把上述问题复杂一下,假设我们现在重复上述实验过程,即第一次,重复实验10次,观察到实验成功次数为1次;第二次,重复实验10次,观察实验成功次数为2次。问:你能告诉我实验成功的次数为几次吗?还是用数学严格的进行求解一次!
这里我们有两个观察值,即随机变量
y1=1,y2=2
,两个随机变量符合相互独立的条件,所以由概率公式得:
同样的,要求 θ 使得似然函数取极大值,我们需要对等式进行求导,问题来了,这是2个观察值,n个观察值进行求导,那这复杂得根本无法计算。因此简单的想法就是把求导的乘法法则能够映射到求导的加法法则,因此便有了对数似然函数的引出,即取 log 函数,得:
这样对上式进行求导便方便很多,更关键的是,求解出来的 θ 值与原先的概率分布函数是等价的。
求得 θ=0.15 。即试验成功的概率为0.15。
6.总结
最大似然估计,只是一种概率论在统计学的应用,它是参数估计的方法之一。说的是已知某个随机样本满足某种概率分布,但是其中具体的参数不清楚,参数估计就是通过若干次实验,观察其结果,利用结果推出参数的大概值。对似然的理解是基于最基本的假设,我们得到的观察结果是在某个参数模型下出现概率最大的。我们并不考虑,实际当实验成功率为0.7时,我们观察两次,分别只观察到10次实验中成功1次和2次。这种情况在某种环境下可能会发生,但倘若发生,我们就认为参数
θ=0.15
,而不是0.7。
求最大似然函数估计值的一般步骤:
1. 写出似然函数
2. 对似然函数取对数,并整理
3. 求导数
4. 解似然方程
上述内容摘自博文–最大似然估计学习总结——MadTurtle,有省略的部分也有补充的部分。
EM算法
1.定义
概率模型有时即含有观测变量,又含有隐变量或潜在变量。如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计法估计模型参数。但是,当模型含有隐变量时,就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法。我们仅讨论极大似然估计,极大后验概率估计与其类似。(隐含了一个问题,概率模型中当存在隐变量时,就无法直接用极大似然估计法进行求解,这是为什么?)
2.三硬币模型
假设有3枚硬币,分别记作A,B,C。这些硬币正面出现的概率分别是
π,p,q
.进行如下投掷试验:先掷硬币A,根据其结果选出硬币B或硬币C,正面选硬币B,反而选硬币C;然后掷选出的硬币,掷硬币的结果,出现正面记作1,出现反面记作0;独立地重复n次试验(这里,n=10),预测结果如下:
假设智能观测到掷硬币的结果,不能观测掷硬币的过程。问如何估计三硬币正面出现的概率,即三硬币模型的参数。
同样的,先用先前似然估计方法来求解一波,看看能否给出答案。假设我们知道了一个观测值:
直接用 θ=(π,p,q) 得
该式子中 y 为已知,其他参数均位置,假设我们知道观察序列的第一次投掷结果为1,因此把
以极大似然方法进行求解,分别对参数 π,p,q 进行求导,你会发现对 π 求导,求出 p=q 来,对 p,q 分别求导求出 π=0和π=1 ,显然是没有解析解。因此,传统的似然估计方法是无法解决上述三硬币模型的问题。这也是为什么EM算法提出的原因,即它能解决传统求导解决不了的问题。(遗留一个问题和一个证明,三硬币模型中是由于 π 的出现,使得似然方程无解?是一旦隐藏变量出现,就无法求解似然方程嘛?如何证明。)
这里,随机变量
y
是观测变量,表示一次试验观测的结果是1或0;随机变量
将观测数据表示为
Y=(Y1,Y2,...,Yn)T
,未观测数据表示为
Z=(Z1,Z2,...,Zn)T
,则观测数据的似然函数为
即
考虑求模型参数 θ=(π,p,q) 的极大似然估计,即
log 函数中有加法,对参数进行求偏导显然是困难的。因此,我们需要另辟蹊径来求解该似然函数的极大值。也就是该似然函数并非是单纯的观测随机变量的概率分布函数,而是隐藏了不可观测变量的概率分布,显然如果并不知道每个样本背后隐藏变量的值,那么求解出来的参数是无意义的。
回到身高分布问题,我们现在假设在全校抽中了100个男生,100个女生,男生和女生的身高都符合正态分布,即背后的概率模型为正态分布概率密度函数,对样本进行统计时,我们标注了他是男生,她是女生,且观测到了每个人的身高,据此我们就可以列出一个关于男生和女生的似然函数,根据似然估计方法,我们通过求导的方式便能求解出模型的参数 θ=(μ,δ)T 。
但假设我们在全校抽中了200个人,但背后我们却没有统计它们的性别,如果把这200个人当作同一群人进行似然函数的建模求解,那显然是不太明智的,因为我们都知道,性别的差异对身高的影响是相当大的,我们非要把女生当男生,并由此计算出模型的参数来,再通过该模型预测出来的男生身高,就显得不那么准确了。所以,在针对隐藏变量的情况下,我们需要考虑性别这样的因素,这也就有了EM算法提出的意义。
在对问题进行分析时,其实我们是可以不用考虑隐藏变量的因素的?但诸如隐马尔可夫模型中的学习问题时,很显然它有隐含的状态,那么进行参数学习时,我们就需要用到EM算法。同样地,三硬币模型也明确的告诉我们含有隐藏变量,因此也必须使用EM算法进行求解。
刚才说了,对含有隐藏变量的似然函数是无法用求导的方式进行求解的,我们先把式子写出来,即
EM算法是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化 L(θ) ,我们可以不断地建立 l 的下界(E步),然后优化下界(M步)。
对于每一个样例
第一步和第二步比较直接,就是分子分母同乘以一个相等的函数。第二步和第三步利用了Jesson不等式,具体的推导过程请参看博文- EM算法原理。
这个过程可以看作是对
L(θ)
求了下界。对于
Qi
的选择,有多种可能,哪种更好呢?假设
θ
给定,那么
L(θ)
的值就决定于
Qi(z(i)
和
p(x(i),z(i);θ)
。我们可以通过调整这两个概率是下界不断上升,以逼近
L(θ)
的真实值,那么什么时候算是调整好了呢?当不等式变成等式时,说明我们调整后的概率都能够等价于
L(θ)
了。按照这个思路,我们要找到等式成立的条件。根据Jessen不等式,要想让等式成立,需要让随机变量变成常数值,这里得到:
C为常数,不依赖于 z(i) 。对此式子做进一步推导,我们知道 ∑ZQi(z(i))=1 ,可以推导得:
对推导过程感兴趣的,可以继续参看博文- EM算法原理。
这里简单提一下取常数的物理含义,上图即为 log 函数,而在 log 函数上两个点,可以分别用 p(x(i),z(i);θ)Qi(z(i)) 和 p(x(j),z(j);θ)Qj(z(j)),i,j 分别表示不同的样本,要让Jessen不等式成立,显然需要让这两点重合,才能取得等号。具体的解释请参看视频教程 七月在线-18分钟理解EM算法。
刚开始接触这个
L(θ)
不理解为什么要让下界函数和
L(θ)
等号成立,然后再开始对新的下界函数进行求极值的过程,并且求极值过程能够逼近
L(θ)
的极值,不急,咱们来看看博文从最大似然到EM算法浅解的解释,参看下图:
E步的过程,就是调整
Q(z)
使得下界
J(Z,Q)
不断上升,直到与
L(θ)
在
θ
点重合,找到了下界函数后,固定
Q(z)
,对参数进行迭代,找到
J(Z,θ)
的极大值,反复上述操作,从而逼近
L(θ)
的极大值。(两个问题,
θ
怎么变?该过程为何收敛?)参看书本《统计学习方法》第159页
算法(EM算法)
输入:观测变量数据Y,隐变量数据Z,联合分布 P(Y,Z|θ) ,条件分布 P(Z|Y,θ)
输出:模型参数 θ
(1) 选择参数的初值 θ(0) ,开始迭代;
(2) E步:记 θ(0) 为第 i 次迭代参数θ 的估计值,在第 i+1 次迭代的E步,计算
Q(θ,θ(i))=∑ZlogP(Y,Z|θ)P(Z|Y,θ(i))
这里, P(Z|Y,θ(i)) 是在给定观测数据Y和当前的参数估计 θ(i) 下隐变量数据Z的条件概率分布。
(3) M步:求使得 Q(θ,θ(i)) 极大化的 θ ,确定第 i+1 次迭代的参数的估计值 θ(i+1)
θ(i+1)=argmaxθQ(θ,θ(i))
(4) 重复第(2)步和第(3)步,直到收敛。
Q函数的推导在书本《统计学习方法》第159页。
行文至此,除了在数学上能够强行“理解”证明的过程,但实际的物理问题却始终无法对应到这些形式化的符号中去,暂且把这些恼人数学放一边,来看看针对某些特定的实际问题,算法是如何解决的,没准能够从中理解一些数学公式的实际含义。
Code Time
双硬币模型
假设有两枚硬币A、B,以相同的概率随机选择一个硬币,进行如下的抛硬币实验:共做5次实验,每次实验独立的抛十次,结果如图中a所示,例如某次实验产生了H、T、T、T、H、H、T、H、T、H,H代表正面朝上。
假设试验数据记录员可能是实习生,业务不一定熟悉,造成a和b两种情况
a表示实习生记录了详细的试验数据,我们可以观测到试验数据中每次选择的是A还是B
b表示实习生忘了记录每次试验选择的是A还是B,我们无法观测实验数据中选择的硬币是哪个
问在两种情况下分别如何估计两个硬币正面出现的概率?
这是实习生a记录的情况,由于这里数据的背后参考模型已知(已分好类),因此用极大似然估计方法就能分别求出
θ^A和θ^B
的概率来。与上文第一节中的例子完全类似。
这是实习生b记录的情况,令人糟糕的是数据的背后参考模型混淆在了一起,我们无法得知它们这几组实验数据是由A抛的还是由B抛的,因为这里隐含了一个该组数据是A还是B的分类问题。抛开数学对隐含变量的复杂求解过程,我们可以先给出一个思路来解决上述问题。
第一,既然我们不清楚是A类还是B类,但假定我们初始化了A类硬币抛正面的概率和B类硬币抛正面的概率,这两者的概率是随意定的,但由于两者概率可以存在差异,假设
P(y=H;θA)>P(y=H;θB)
,那么一个明显的特征就是,由于能观察到10次硬币中有几次是成功的,我们可以基于这次观察,给出
P(z=A|Y;θA,θB)
的概率,上式的含义是可以根据两个参数的初值求出,在给定观察序列的情况下,它属于A类还是B类的概率。用公式可以表示为:
其中,z表示单个观察到的随机变量,此处z=A or B(属于分类问题),Y表示观察序列,即 Y=(y1,y2,...,y10)T 。由此,给定观察序列后,我们可以算出属于A类的概率和属于B类的概率,那么很显然CoinA 和CoinB 不再是属于你有我没有,你没有我有的敌对关系,因为我自己都还不是很清楚是不是A类,由此10个硬币,我们就根据概率进行一次平均分配咯,这样CoinA 和CoinB 在一次观察结果时,都能得到属于自己的那一份,非常的和谐。这一部分便是求期望的过程,即对于第一个观察序列中,10次抛硬币过程中5次为正面朝上,令 yj=5 ,由此可以得到关于隐含变量的数学期望 E(z)=0.45∗5+0.55∗5 ,“+”号的左边右边,分别表示CoinA的分配和CoinB的分配。分配的份额根据z函数的分布给定,z函数的分布规则根据缺失的信息进行建模,解由初始参数求出。
因此分类问题,给出了每个CoinA 和CoinB 的分配额,有了所有观察值CoinA和CoinB的分配额,我们就可以单独的对CoinA和CoinB进行最大似然估计方法。求出来的新参数,再求z函数,求期望,求参数,如此迭代下去,直到收敛到某一结果。
算法实现
在双硬币模型中,对某个种类硬币投掷10次中成功n次概率模型 P(y|θ)=(10n)θn(1−θ)(10−n),y=0,1,...,10 符合伯努利分布。
前期准备工作
我们可以实际的操作一把,来看看在python中如何可视化伯努利分布。前期环境搭建:
1. 安装python扩展包:scipy,matplotlib,seaborn,ipywidgets。
- 在安装扩展包过程中遇到了一些课,如安装scipy,用pip直接安装遇到安装失败的情况,如果有类似的情况可以参考链接[Python]Windows7 x64安装numpy和scipy
2. 安装ipython,貌似在python3.5版本中,自带了ipython。
- 安装完毕,在cmd中敲入ipython notebook,开始编程之旅
进入游览器
点击new按钮后,新增python3按钮,进入下一界面
实际编写
接下类就可以进行我们实际的操作了,首先引入一些必须的科学计算包:
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
from ipywidgets import interact, FloatSlider
%matplotlib inline
编写伯努利分布函数
a=range(11)
def plot_binomial(p=0.5):
fig, ax = plt.subplots(figsize=(4,3))
y = [0]*11
for i in a:
y[i-1] = stats.binom.pmf(i, 10, p)
ax.bar(a,y,label="$p = %.1f$" % p)
ax.set_ylabel('PMF of $k$ heads')
ax.set_xlabel('$k$')
ax.set_ylim((0,0.5))
ax.set_xlim((0,10))
ax.legend()
return fig
可视化伯努利分布函数
interact(plot_binomial, p=FloatSlider(min=0.0,max=1.0,step=0.1,value=0.2))
可视化界面如下
这是在p=0.2的情况下的伯努利分布函数,代回双硬币模型中去,当观察到10次实验中只有2次成功了,那么该
θ
参数便是0.2。因为只有当
θ=0.2
时,10次实验中出现成功次数为2次的概率最大。
回到实习生b,由数据可得观察矩阵为:
observations = np.array([[1,0,0,0,1,1,0,1,0,1],
[1,1,1,1,0,1,1,1,1,1],
[1,0,1,1,1,1,1,0,1,1],
[1,0,1,0,0,0,1,1,0,0],
[0,1,1,1,0,1,1,1,0,1]])
实际每组观察数据属于A类,B类的隐藏状态为:
coins_id = np.array([False,True,True,False,True])
那么在观察数组中,属于A类的数组为:
observations[coins_id]
输出结果为:
array([[1, 1, 1, 1, 0, 1, 1, 1, 1, 1],
[1, 0, 1, 1, 1, 1, 1, 0, 1, 1],
[0, 1, 1, 1, 0, 1, 1, 1, 0, 1]])
在所有属于A类的数组中,总的实验成功次数为:
np.sum(observations[coins_id])
输出结果为:
24
所以说,针对属于A类的硬币,它的参数 θA :
1.0*np.sum(observations[coins_id])/observations[coins_id].size
输出结果为:
0.80000000000000004
同理,对于属于B类的硬币,它的参数为 θB :
1.0*np.sum(observations[~coins_id])/observations[~coins_id].size
输出结果为:
0.45000000000000001
EM算法步骤
首先来看看,针对第一组观察数据,每一步的数据是如何求出的。
# 对于实验成功率为0.6的情况,10次中成功5次的概率
coin_A_pmf_observation_1= stats.binom.pmf(5,10,0.6)
coin_A_pmf_observation_1
# 均是针对第一组观察数据的情况
coin_B_pmf_observation_1= stats.binom.pmf(5,10,0.5)
coin_B_pmf_observation_1
# 针对第一组观察数据中,属于A类硬币投掷的概率p(z=A|Y,theta)
normalized_coin_A_pmf_observation_1 = coin_A_pmf_observation_1/(coin_A_pmf_observation_1+coin_B_pmf_observation_1)
print ("%0.2f" %normalized_coin_A_pmf_observation_1)
# 针对第一组观察数据中,属于B类硬币投掷的概率p(z=B|Y,theta)
normalized_coin_B_pmf_observation_1 = coin_B_pmf_observation_1/(coin_A_pmf_observation_1+coin_B_pmf_observation_1)
print ("%0.2f" %normalized_coin_B_pmf_observation_1)
#求期望过程
weighted_heads_A_obervation_1 = 5*normalized_coin_A_pmf_observation_1
print ("Coin A Weighted count for heads in observation 1: %0.2f" %weighted_heads_A_obervation_1)
weighted_tails_A_obervation_1 = 5*(1-normalized_coin_A_pmf_observation_1)
print ("Coin A Weighted count for tails in observation 1: %0.2f" %weighted_tails_A_obervation_1)
weighted_heads_B_obervation_1 = 5*normalized_coin_B_pmf_observation_1
print ("Coin B Weighted count for heads in observation 1: %0.2f" %weighted_heads_B_obervation_1)
weighted_tails_B_obervation_1 = 5*(1-normalized_coin_B_pmf_observation_1)
print ("Coin B Weighted count for tails in observation 1: %0.2f" %weighted_tails_B_obervation_1)
单步EM算法,迭代一次算法实现步骤
def em_single(priors,observations):
"""
performs a single Em step
Arguments
---------
priors:[theta_A,theta_B]
observations:[m X n matrix]
Returns
-------
new_priors:[new_theta_A,new_theta_B]
"""
counts ={'A':{'H':0,'T':0},'B':{'H':0,'T':0}}
theta_A = priors[0]
theta_B = priors[1]
# E setp
for observation in observations:
len_observation = len(observation)
# 对数组求和 head =1,tail =0
num_heads = observation.sum()
num_tails = len_observation - num_heads
# 如果属于A类,A的分布概率 和属于B类的分布概率
contribution_A = stats.binom.pmf(num_heads,len_observation,theta_A)
contribution_B = stats.binom.pmf(num_heads,len_observation,theta_B)
# z分布概率
weight_A = contribution_A /(contribution_A+contribution_B)
weight_B = contribution_B /(contribution_A+contribution_B)
# Incrementing counts
counts['A']['H'] += weight_A*num_heads
counts['A']['T'] += weight_A*num_tails
counts['B']['H'] += weight_B*num_heads
counts['B']['T'] += weight_B*num_tails
# M setp
new_theta_A = counts['A']['H']/(counts['A']['H']+counts['A']['T'])
new_theta_B = counts['B']['H']/(counts['B']['H']+counts['B']['T'])
return [new_theta_A, new_theta_B]
迭代一次输出结果为:
em_single([0.6,0.5],observations)
# 结果
[0.71301223540051617, 0.58133930831366265]
收敛算法
def em(observations,prior,tol=1e-6,iterations=1000):
import math
iteration = 0
while iteration <iterations:
new_prior = em_single(prior,observations)
delta_change = np.abs(prior[0] -new_prior[0])
if delta_change < tol:
break
else:
prior = new_prior
iteration +=1
return [new_prior,iteration]
最终结果为:
em(observations, [0.6,0.5])
# 结果
[[0.79678875938310978, 0.51958393567528027], 14]
详细代码请参看博文Programatically understanding Expectation Maximization algorithm
最终实习生b的EM算法得出的结果,跟实习生a得出的参数还是非常相近的,但EM算法跟初始值的设置有着很大的关系,不信,修改[06,0.5]多尝试尝试。
针对EM算法的所有相关内容都已经阐述完毕了,在学习过程中,遇到了许多坑,诸如数学公式与实际的问题没法完全一一对应,但个人认为,自家数学功底不够深厚时,唯有参照实际情况,来解决一个实际问题,多加练习,在日后回过头来看数学定义时,没准能够恍然大悟。解决双硬币模型的思想是很简单的,但对当细节问题进行深入理解时,本文还有很多不足,如算法如何收敛,为何收敛,算法E-step为何是求期望。待补足数学上的漏洞或者经历过大量实践后,有望能够解决,未完待续。