0、前言
EM算法,也叫最大期望算法,或者是期望最大化算法,是机器学习十大算法之一,它很简单,但是也同样很有深度。
-
简单是因为它就分两步求解问题:
- E步:求期望(expectation)
- M步:求极大(maximization)
-
深度在于它的数学推理涉及到比较繁杂的概率公式等。
所以本文会介绍很多概率方面的知识,不懂的同学可以先去了解一些知识,当然本文也会尽可能的讲解清楚这些知识,讲的不好的地方麻烦大家评论指出,后续不断改进完善。
1、EM算法引入
概率模型有时候既含有观测变量,又含有隐变量或潜在变量,如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计方法估计模型参数,但是,当模型含有隐变量时,就不能简单的使用这些方法,EM算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法,我们讨论极大似然估计,极大后验概率估计与其类似。
参考统计学习方法书中的一个例子来引入EM算法:
假设有3枚硬币,分别记做A、B、C,这些硬币正面出现的概率分别是
π、
p、
q,进行如下实验:
- 先掷硬币A,根据硬币A的结果选出硬币B和硬币C,正面选硬币B,反面选硬币C;
- 通过选择出的硬币(B或者C),记录掷硬币的结果出现正面为1,反面为0。
也就是说一共需要掷两次硬币,第一次为硬币A,第二次硬币需要根据A的结果来选择,只记录一次结果,就是第二次掷的硬币的正反面。把这整个过程看作一次实验,如此独立地重复n次实验。我们当前规定n=10,则10次的结果如下所示:
1,1,0,1,0,0,1,0,1,1
假设现在我们只知道到记录掷硬币的结果,也就是上面的10次结果,不知道每次实验掷硬币的过程,也就是不知道每次实验第一次掷硬币A的结果;现在需要估计每次实验结果出现正面的概率,怎么求呢?
我们来构建这样一个三硬币模型:
硬币A的结果 |
第二次的硬币 |
第二次的结果 |
概率 |
正面 |
硬币B |
正面 |
πp |
正面 |
硬币B |
反面 |
π(1−p) |
反面 |
硬币C |
正面 |
πq |
反面 |
硬币C |
反面 |
π(1−q) |
记每次实验记录的结果为
y,则第二次硬币为B的概率为
π,则B的结果概率为
πpy(1−p)(1−y),当出现正面,就令
y=1,概率就为
πp。对C同理。计算模型的概率
P(y∣θ),其中
θ=(π,p,q)是模型的参数,
y受到硬币A的结果的影响,记硬币A的结果为
z,则
P(y∣θ)是联合概率
P(y,z∣θ)的边际概率,即
P(y∣θ)=∑zP(y,z∣θ),所以:
P(y∣θ)=z∑P(y,z∣θ)=z∑P(z∣θ)P(y∣z,θ)=P(z=正面∣θ)P(y∣z=正面,θ)+P(z=反面∣θ)P(y∣z=反面,θ)=πpy(1−p)(1−y)+(1−π)qy(1−q)(1−y)
- 若
y=1,表示这此看到的是正面,这个正面有可能是B的正面,也可能是C的正面,则
P(y=1∣θ)=πp+(1−π)q
- 若
y=0,则
P(0∣θ)=π(1−p)+(1−π)(1−q)
y是观测变量,表示一次观测结果是
1或
0,
z是隐藏变量,表示掷硬币A的结果,这个是观测不到结果的,
θ=(π,p,q)表示模型参数。若考虑
n次实验,将观测数据表示为
Y=(Y1,Y2,...,Yn)T,未观测的数据表示为
Z=(Z1,Z2,...,Zn)T,则观测函数的似然函数是:
P(Y∣θ)=Z∑P(Z∣θ)P(Y∣Z,θ)=i=0∏n[πpyi(1−p)(1−yi)+(1−π)qyi(1−q)(1−yi)]
考虑求模型参数
θ=(π,p,q)的极大似然估计,即:
θ^=argθmaxlogP(Y∣θ)
这个问题没有解析解,只有通过迭代方法来求解,EM算法就是可以用于求解这个问题的一种迭代算法,下面给出EM算法的迭代过程:
-
首先选取初始值,记做
θ0=(π0,p0,q0),第
i次的迭代参数的估计值为
θi=(πi,pi,qi)
-
E步:计算在模型参数
πi,pi,qi下观测变量
yi来源于硬币B的概率:
P(zi=πi∣yi,θi)=μi+1=P(yi∣θi)P(zi=πi,yi∣θi)=zi∑P(zi,yi∣θi)P(zi=πi,yi∣θi)=P(zi=πi,yi∣θi)+P(zi=(1−πi),yi∣θi)P(zi=πi,yi∣θi)=πi(pi)yi(1−pi)1−yi+(1−πi)(qi)yi(1−pi)1−yiπi(pi)yi(1−pi)1−yi
备注一下:这个公式的分母是
P(Y∣θ),分子表示结果来源于B硬币的概率
P(Y,Z∣θ)。
-
M步:计算模型参数的新估计值:
πi+1=n1j=1∑nμji+1
因为第二次使用B硬币是基于A硬币出现正面的结果,所以A硬币出现正面的概率就是
μj的平均值。
pi+1=j=1∑nμji+1j=1∑nμji+1yj
分子乘以
yi,所以其实是计算B硬币出现正面的概率。
qi+1=j=1∑n(1−μji+1)j=1∑n(1−μji+1)yj
其中,
(1−μji+1)表示出现C硬币的概率,即A出现反面的概率。
闭环形成,从
P(Y∣θ) 到
π、p、q一个闭环流程,接下来可以通过迭代法来做完成。针对上述例子,我们假设初始值为
π0=0.5,p0=0.5,q0=0.5,因为对
yi=1和
yi=0均有
μj1=0.5,利用迭代公式计算得到
π1=0.5,p1=0.6,q1=0.6,继续迭代得到最终的参数:
π0
=0.5,p0
=0.6,q0
=0.6
如果一开始初始值选择为:
π0=0.4,p0=0.6,q0=0.7,那么得到的模型参数的极大似然估计是:
π
=0.4064,p
=0.5368,q
=0.6432
这说明EM算法与初值的选择有关,选择不同的初值可能得到不同的参数估计值。
这个例子中你只观察到了硬币抛完的结果,并不了解A硬币抛完之后,是选择了B硬币抛还是C硬币抛,这时候概率模型就存在着隐含变量!
2、具体的EM算法
-
输入:观测变量数据Y,隐变量数据Z,联合分布
P(Y,Z∣θ),条件分布
P(Z∣Y,θ);
-
输出:模型参数
θ;
-
(1) 选择参数的初值
θ0,开始迭代
-
(2) E步:记
θi为第
i次迭代参数
θ的估计值,在第
i+1次迭代的E步,计算
Q(θ,θi)函数:
Q(θ,θi)=EZ[logP(Y,Z∣θ)∣Y,θi]=Z∑logP(Y,Z∣θ)P(Z∣Y,θi)
这里,
P(Z∣Y,θi)是在给定观测数据
Y和当前的参数估计
θi下隐变量数据
Z的条件概率分布;
- (3) M步:求使
Q(θ,θi)极大化的
θ,确定第i+1次迭代的参数的估计值
θi+1,
θi+1=argθmaxQ(θ,θi)
其中,
Q(θ,θi)是EM算法的核心,称为Q函数(Q function),这个是需要自己构造的,可以看出这是一个期望,而M步是极大化这个期望值,这也是EM(Expectation Maximization )算法名字的来源。
- (4) 重复第(2)步和第(3)步,直到收敛,收敛条件:
∣∣θi+1−θi∣∣<ε1或者:∣∣Q(θi+1,θi)−Q(θi,θi)∣∣<ε2
收敛迭代就结束了。我们来拆解一下这个M步骤,
3、EM算法推导
3.1 Jensen不等式
Jensen不等式,这个公式在算法推导和证明收敛都能用到,Jensen不等式主要是如下的结论:
f(E(X))≤E(f(x))
f(E(X))≥E(f(x))
3.2 EM推导
如果模型不包含隐变量,我们在模型学习的时候,直接使用极大似然估计使
P(Y∣θ)最大,其中Y是观测数据,
θ是模型参数,意思就是在参数
θ的情况下,观测数据Y出现的概率最大。然而当我们面对一个含有隐变量的概率模型时,我们也想使在参数
θ的情况下,观测数据Y出现的概率最大,即最大化
P(Y∣θ)。但是我们不能直接使用极大似然估计,因为参数
θ包含有隐变量的信息,而Y与隐变量相关,估计Y需要知道Z,估计Z需要知道Y,形成一个死循环。可以将隐变量分离出来,对于每一个样本,我们需要确定Z来使联合概率
P(Y,Z∣θ)最大,由条件概率公式:
P(Y,Z∣θ)=P(Y∣Z,θ)P(Z∣θ):
P(Y∣θ)=Z∑P(Y,Z∣θ)=Z∑P(Y∣Z,θ)P(Z∣θ)
然后就可以使用极大似然估计使
P(Y∣θ)出现的概率最大,为了方便计算,对似然函数取对数。所以我们的目标是极大化观测数据Y关于参数
θ的对数似然函数
L(θ):
L(θ)=logP(Y∣θ)=logZ∑P(Y∣Z,θ)P(Z∣θ)
想要直接极大化上面的式子,可以对Y,Z求偏导,但是log里面连加求导,会很麻烦,可以自己想象。可以看Jensen不等式正好将求和和函数连在一起,但是必须要使求和的那个系数之和为1,刚好和概率之和为1对上了呀。构造一个概率之和,而
Z∑P(Z∣Y,θ)=1,刚好可以利用Jensen不等式。
根据EM算法的定义,我们需要通过迭代的方式逐渐去极大化
L(θ),每一次更新
θ时,都能使
L(θ)增大,也就是假设当前第
i次迭代的参数估计值为
θi,则在这
i+1次迭代中,我们需要重新估计参数值,使
L(θ)增大,即
L(θ)>L(θi)。则
L(θ)−L(θi):
L(θ)−L(θi)(构造Jensen不等式←)(Z∑P(Z∣Y,θi)=1←)(log是凸函数←)=logZ∑P(Y∣Z,θ)P(Z∣θ)−logP(Y∣θi)=log(Z∑P(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi)P(Z∣Y,θi))−logP(Y∣θi)=log(Z∑P(Z∣Y,θi)P(Z∣Y,θi)P(Y∣Zθ)P(Z∣θ))−logP(Y∣θi)≥Z∑P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θi)
这一个步骤就是采用了凸函数的Jensen不等式做转换。因为
Z是隐藏变量,所以有:
Z∑P(Z∣Y,θi)=1,P(Z∣Y,θi)>0
于是继续变:
L(θ)−L(θi)=log(Z∑P(Z∣Y,θi)P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θi)≥Z∑P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θi)=Z∑P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ))−Z∑P(Z∣Y,θi)logP(Y∣θi)=Z∑P(Z∣Y,θi)[logP(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)−logP(Y∣θi)]=Z∑P(Z∣Y,θi)logP(Z∣Y,θi)P(Y∣θi)P(Y∣Z,θ)P(Z∣θ)
也就是:
L(θ)≥L(θi)+Z∑P(Z∣Y,θi)log(P(Y∣Z,θi)P(Y∣θi)P(Y∣Z,θ)P(Z,θ)),有下界(不等式右边就是下界),任何使下界增大的
θ都可以使
L(θ)增大,因此我们需要最大化下界,来得到近似的极大值。
第
i+1次迭代的参数估计值
θi+1,就是使下界最大化的
θ:
θi+1=L(θi)+Z∑P(Z∣Y,θi)log(P(Y∣Z,θi)P(Y∣θi)P(Y∣Z,θ)P(Z∣θ))=L(θi)+Z∑P(Z∣Y,θi)logP(Y∣Z,θi)P(Y∣Z,θ)P(Z∣θ)−L(θi)=argθmaxZ∑P(Z∣Y,θi)logP(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)=argθmaxZ∑P(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)−Z∑P(Z∣Y,θi)logP(Z∣Y,θi))=argθmaxZ∑P(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)=argθmaxZ∑P(Z∣Y,θi)logP(Y,Z∣θ)=argθmaxQ(θ,θi)
其中
log分母提出来是关于
Z的
Z∑P(Z∣Y,θi)logP(Z∣Y,θi),是一个常数可以去掉。其中
Q(θ,θi)=Z∑P(Z∣Y,θi)logP(Y,Z∣θ)是EM算法中非常重要的Q函数,很多时候写成期望的形式:
Q(θ,θi)=Z∑P(Z∣Y,θi)logP(Y,Z∣θ)=EZ[logP(Y,Z∣θ)∣Y,θi]=EZ∣Y,θilogP(Y,Z∣θ)
记
B(θ,θi)=L(θi)+Z∑P(Z∣Y,θi)log(P(Y∣Z,θi)P(Y∣θi)P(Y∣Z,θ)P(Z,θ)),也就是
L(θ)的下界,可以看出
L(θ)≥B(θ,θi),当
θ=θi时,等号成立。函数
B(θ,θi)的增加会导致
L(θ)在每次迭代都是增加的,我们想尽快迭代到最优值,就要求在每一次迭代的时候都取到函数
B(θ,θi)的最大值,也就是
θi+1,也就是Q函数的最大值。然后基于
θi+1,重新计算
B(θ,θi+1),会得到类似的图,继续寻找最大值
θi+2,以此类推下去。在这个过程中,对数似然函数
L(θ)不断增大,直到达到全局最优值。
3.3 EM算法的收敛性
EM算法的最大优点是简单性和普适性,但是EM算法得到的估计序列是否收敛呢?是收敛到全局最大值还是局部最大值?
定理一:设
P(Y∣θ)是观测数据的似然函数,
θi(i=1,2,...)是EM算法得到的参数估计序列,
P(Y∣θi)(i=1,2,...)是对应的似然函数序列,则
P(Y∣θi)是单调递增的,即:
P(Y∣θi+1)≥P(Y∣θi)
证明过程:
由条件概率
P(Y,Z∣θ)=P(Y∣θ)P(Z∣Y,θ)可得:
P(Y∣θ)=P(Z∣Y,θ)P(Y,Z∣θ)
取对数有:
logP(Y∣θ)=logP(Y,Z∣θ)−logP(Z∣Y,θ)
由上文可知:
Q(θ,θi)=Z∑P(Z∣Y,θi)logP(Y,Z∣θ)
令:
H(θ,θi)=Z∑P(Z∣Y,θi)logP(Z∣Y,θ)
Q(θ,θi)−H(θ,θi)=Z∑P(Z∣Y,θi)logP(Y,Z∣θ)−Z∑P(Z∣Y,θi)logP(Z∣Y,θ)=Z∑P(Z∣Y,θi)logP(Z∣Y,θ)P(Y,Z∣θ)=Z∑P(Z∣Y,θi)logP(Y∣θ)=logP(Y∣θ)
则:
logP(Y∣θi+1)−logP(Y∣θi)=Q(θi+1,θi)−H(θi+1,θi)−[Q(θi,θi)−H(θi,θi)]=Q(θi+1,θi)−Q(θi,θi)−[H(θi+1,θi)−H(θi,θi)]
要证明
P(Y∣θi+1)≥P(Y∣θi),只需证明上式的右边是非负的。其中
Q(θi+1,θi)−Q(θi,θi)是必大于0的,因为
θi+1是
Q(θ,θi)的最大值。则:
H(θi+1,θi)−H(θi,θi)(利用Jensen不等式←)=Z∑P(Z∣Y,θi+1)logP(Z∣Y,θ)−Z∑P(Z∣Y,θi)logP(Z∣Y,θ)=Z∑P(Z∣Y,θi)logP(Z∣Y,θi)P(Z∣Y,θi+1)=log(Z∑P(Z∣Y,θi)P(Z∣Y,θi)P(Z∣Y,θi+1))=log(Z∑P(Z∣Y,θi+1))=log1=0
由此可以证明
P(Y∣θi+1)≥P(Y∣θi)。
定理二:设
L(θ)=logP(Y∣θ)为观测数据的对数似然函数,
θi(i=1,2,...)是EM算法得到的参数估计序列,
L(θi)(i=1,2,...)是对应的对数似然函数序列。
- 如果
P(Y∣θ)有上界,则$
L(θi)=logP(Y∣θi)收敛到某一值
L∗;
- 在函数
Q(θ,θ′)与
L(θ)满足一定条件下,由EM算法得到的参数估计序列
θi的收敛值
θ∗是
L(θ)的稳定点。
这里就不证明了
4、EM算法在高斯混合模型中的应用
4.1 高斯混合模型
EM算法的一个重要应用场景就是高斯混合模型的参数估计。高斯混合模型就是由多个高斯模型组合在一起的混合模型(可以理解为多个高斯分布函数的线性组合,理论上高斯混合模型是可以拟合任意类型的分布),例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的:
两个高斯模型可以拟合数据集,如图所示:
如果有多个高斯模型,公式表示为:
P(y∣θ)=k=1∑Kakϕ(y∣θk)ϕ(y∣θk)=2π
δk1exp(−2δk2(y−μk)2)ak>0,∑ak=1
ϕ(y∣θk)表示为第k个高斯分布密度模型,定义如上,其中
ak表示被选中的概率。在本次模型
P(y∣θ)中,观测数据是已知的,而观测数据具体来自哪个模型是未知的,有点像之前提过的三硬币模型,我们来对比一下,A硬币就像是概率
ak,用来表明具体的模型,而B、C硬币就是具体的模型,只不过这里有很多个模型,不仅仅是B、C这两个模型。我们用
γjk来表示,则:
γjk={10第j个观测数据来源于第k个模型否则
所以一个观测数据
yj的隐藏数据
(γj1,γj2,...,γjk),那么完全似然函数就是:
P(y,γ∣θ)=k=1∏Kj=1∏N[akϕ(y∣θk)]γjk
取对数之后等于:
log(P(y,γ∣θ))=log(k=1∏Kj=1∏N[akϕ(y∣θk)]γjk)=K∑k=1(j=1∑N(γjk)log(ak)+j=1∑N(γjk)[log(2π
1)−log(δk)−2δk2(yi−μk)2])
- E 步 :
Q(θ.θi)=E[log(P(y,γ∣θ))]=K∑k=1(j=1∑N(Eγjk)log(ak)+j=1∑N(Eγjk)[log(2π
1)−log(δk)−2δk2(yi−μk)2])
其中我们定义
γjk^:
γjk^=E(γjk∣y,θ)=∑k=1Kakϕ(yi∣θk)akϕ(yi∣θk)j=1,2,..,N;k=1,2,...,Knk=j=i∑NEγjk
于是化简得到:
Q(θ.θi)=K∑k=1(nklog(ak)+j=1∑N(Eγjk)[log(2π
1)−log(δk)−2δk2(yi−μk)2])
E 步 在代码设计上只有
γjk^有用,用于M步的计算。
- M步,
θi+1=argθmaxQ(θ,θi)
对
Q(θ,θi)求导,得到每个未知量的偏导,使其偏导等于0,求解得到:
μk^=∑j=1Nγjk^∑j=1Nγjk^yiδk^=∑j=1Nγjk^∑j=1Nγjk^(yi−μk)2ak^=N∑j=1Nγjk^
给一个初始值,来回迭代就可以求得值内容。这一块主要用到了
Q(θ.θi)的导数,并且用到了E步的
γjk^。
4.2 混合高斯分布模型python实现
EM算法更多是一种思想,用概率来解决问题的一种方法,具体的代码看自己选用模型,所以并没有通用的模型,本此代码主要是讲解混合高斯分布模型的
这其中的M步 完全按照了 公式来计算。
import numpy as np
import random
import math
import time
'''
数据集:伪造数据集(两个高斯分布混合)
数据集长度:1000
------------------------------
运行结果:
----------------------------
the Parameters set is:
alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0
----------------------------
the Parameters predict is:
alpha0:0.4, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9
----------------------------
'''
def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):
'''
初始化数据集
这里通过服从高斯分布的随机函数来伪造数据集
:param mu0: 高斯0的均值
:param sigma0: 高斯0的方差
:param mu1: 高斯1的均值
:param sigma1: 高斯1的方差
:param alpha0: 高斯0的系数
:param alpha1: 高斯1的系数
:return: 混合了两个高斯分布的数据
'''
length = 1000
data0 = np.random.normal(mu0, sigma0, int(length * alpha0))
data1 = np.random.normal(mu1, sigma1, int(length * alpha1))
dataSet = []
dataSet.extend(data0)
dataSet.extend(data1)
random.shuffle(dataSet)
return dataSet
def calcGauss(dataSetArr, mu, sigmod):
'''
根据高斯密度函数计算值
依据:“9.3.1 高斯混合模型” 式9.25
注:在公式中y是一个实数,但是在EM算法中(见算法9.2的E步),需要对每个j
都求一次yjk,在本实例中有1000个可观测数据,因此需要计算1000次。考虑到
在E步时进行1000次高斯计算,程序上比较不简洁,因此这里的y是向量,在numpy
的exp中如果exp内部值为向量,则对向量中每个值进行exp,输出仍是向量的形式。
所以使用向量的形式1次计算即可将所有计算结果得出,程序上较为简洁
:param dataSetArr: 可观测数据集
:param mu: 均值
:param sigmod: 方差
:return: 整个可观测数据集的高斯分布密度(向量形式)
'''
result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))
return result
def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):
'''
EM算法中的E步
依据当前模型参数,计算分模型k对观数据y的响应度
:param dataSetArr: 可观测数据y
:param alpha0: 高斯模型0的系数
:param mu0: 高斯模型0的均值
:param sigmod0: 高斯模型0的方差
:param alpha1: 高斯模型1的系数
:param mu1: 高斯模型1的均值
:param sigmod1: 高斯模型1的方差
:return: 两个模型各自的响应度
'''
gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)
sum = gamma0 + gamma1
gamma0 = gamma0 / sum
gamma1 = gamma1 / sum
return gamma0, gamma1
def M_step(muo, mu1, gamma0, gamma1, dataSetArr):
mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)
mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)
sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0))
sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1))
alpha0_new = np.sum(gamma0) / len(gamma0)
alpha1_new = np.sum(gamma1) / len(gamma1)
return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new
def EM_Train(dataSetList, iter=500):
'''
根据EM算法进行参数估计
算法依据“9.3.2 高斯混合模型参数估计的EM算法” 算法9.2
:param dataSetList:数据集(可观测数据)
:param iter: 迭代次数
:return: 估计的参数
'''
dataSetArr = np.array(dataSetList)
alpha0 = 0.5
mu0 = 0
sigmod0 = 1
alpha1 = 0.5
mu1 = 1
sigmod1 = 1
step = 0
while (step < iter):
step += 1
gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)
return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1
if __name__ == '__main__':
start = time.time()
alpha0 = 0.3
mu0 = -2
sigmod0 = 0.5
alpha1 = 0.7
mu1 = 0.5
sigmod1 = 1
dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)
print('---------------------------')
print('the Parameters set is:')
print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
))
alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)
print('----------------------------')
print('the Parameters predict is:')
print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
))
print('----------------------------')
print('time span:', time.time() - start)
E步主要计算内容:
其中我们定义
γjk^:
γjk^=E(γjk∣y,θ)=∑k=1Kakϕ(yi∣θk)akϕ(yi∣θk)j=1,2,..,N;k=1,2,...,Knk=j=i∑NEγjk
M步 主要计算内容
这一步骤主要是对Q函数求导后的数据进行计算,利用了 E 步 的
γjk^
import math
import copy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def generate_data(sigma,N,mu1,mu2,mu3,mu4,alpha):
global X
X = np.zeros((N, 2))
X=np.matrix(X)
global mu
mu = np.random.random((4,2))
mu=np.matrix(mu)
global excep
excep=np.zeros((N,4))
global alpha_
alpha_=[0.25,0.25,0.25,0.25]
for i in range(N):
if np.random.random(1) < 0.1:
X[i,:] = np.random.multivariate_normal(mu1, sigma, 1)
elif 0.1 <= np.random.random(1) < 0.3:
X[i,:] = np.random.multivariate_normal(mu2, sigma, 1)
elif 0.3 <= np.random.random(1) < 0.6:
X[i,:] = np.random.multivariate_normal(mu3, sigma, 1)
else:
X[i,:] = np.random.multivariate_normal(mu4, sigma, 1)
print("可观测数据:\n",X)
print("初始化的mu1,mu2,mu3,mu4:",mu)
def e_step(sigma,k,N):
global X
global mu
global excep
global alpha_
for i in range(N):
denom=0
for j in range(0,k):
denom += alpha_[j]* math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:])) /np.sqrt(np.linalg.det(sigma))
for j in range(0,k):
numer = math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/np.sqrt(np.linalg.det(sigma))
excep[i,j]=alpha_[j]*numer/denom
print("隐藏变量:\n",excep)
def m_step(k,N):
global excep
global X
global alpha_
for j in range(0,k):
denom=0
numer=0
for i in range(N):
numer += excep[i,j]*X[i,:]
denom += excep[i,j]
mu[j,:] = numer/denom
alpha_[j]=denom/N
def plotShow():
plt.subplot(221)
plt.scatter(X[:,0].tolist(), X[:,1].tolist(),c='b',s=25,alpha=0.4,marker='o')
plt.title('random generated data')
plt.subplot(222)
plt.title('classified data through EM')
order=np.zeros(N)
color=['b','r','k','y']
for i in range(N):
for j in range(k):
if excep[i,j]==max(excep[i,:]):
order[i]=j
probility[i] += alpha_[int(order[i])]*math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/(np.sqrt(np.linalg.det(sigma))*2*np.pi)
plt.scatter(X[i, 0], X[i, 1], c=color[int(order[i])], s=25, alpha=0.4, marker='o')
ax = plt.subplot(223, projection='3d')
plt.title('3d view')
for i in range(N):
ax.scatter(X[i, 0], X[i, 1], probility[i], c=color[int(order[i])])
plt.show()
if __name__ == '__main__':
iter_num=1000
N=500
k=4
probility = np.zeros(N)
u1=[5,35]
u2=[30,40]
u3=[20,20]
u4=[45,15]
sigma=np.matrix([[30, 0], [0, 30]])
alpha=[0.1,0.2,0.3,0.4]
generate_data(sigma,N,u1,u2,u3,u4,alpha)
for i in range(iter_num):
err=0
err_alpha=0
Old_mu = copy.deepcopy(mu)
Old_alpha = copy.deepcopy(alpha_)
e_step(sigma,k,N)
m_step(k,N)
print("迭代次数:",i+1)
print("估计的均值:",mu)
print("估计的混合项系数:",alpha_)
for z in range(k):
err += (abs(Old_mu[z,0]-mu[z,0])+abs(Old_mu[z,1]-mu[z,1]))
err_alpha += abs(Old_alpha[z]-alpha_[z])
if (err<=0.001) and (err_alpha<0.001):
print(err,err_alpha)
break
plotShow()
5 EM算法在HMM模型中的应用
5.1 HMM模型
HMM模型也是一个含有隐变量的模型,具体参考HMM模型详解
5.2 EM算法在HMM中应用
HMM模型在学习参数的时候,可以使用极大似然估计,也可以使用EM算法(在HMM称为bamu-welch算法)
class HMM(object):
def __init__(self, n, m, a=None, b=None, pi=None):
self.N = n
self.M = m
self.A = a
self.B = b
self.Pi = pi
self.X = None
self.Y = None
self.T = 0
self.alpha = None
self.beta = None
def baum_welch(self, x, criterion=0.001):
self.T = len(x)
self.X = x
while True:
_ = self.forward(self.X)
_ = self.backward(self.X)
xi = np.zeros((self.T - 1, self.N, self.N), dtype=float)
for t in range(self.T - 1):
denominator = np.sum(np.dot(self.alpha[t, :], self.A) *
self.B[:, self.X[t + 1]] * self.beta[t + 1, :])
for i in range(self.N):
molecular = self.alpha[t, i] * self.A[i, :] * self.B[:, self.X[t+1]]*self.beta[t+1, :]
xi[t, i, :] = molecular / denominator
gamma = np.sum(xi, axis=2)
prod = (self.alpha[self.T-1, :]*self.beta[self.T-1, :])
gamma = np.vstack((gamma, prod / np.sum(prod)))
new_pi = gamma[0, :]
new_a = np.sum(xi, axis=0) / np.sum(gamma[:-1, :], axis=0).reshape(-1, 1)
new_b = np.zeros(self.B.shape, dtype=float)
for k in range(self.B.shape[1]):
mask = self.X == k
new_b[:, k] = np.sum(gamma[mask, :], axis=0) / np.sum(gamma, axis=0)
if np.max(abs(self.Pi - new_pi)) < criterion and \
np.max(abs(self.A - new_a)) < criterion and \
np.max(abs(self.B - new_b)) < criterion:
break
self.A, self.B, self.Pi = new_a, new_b, new_pi
参考
主要参考统计学习方法
统计学习方法-代码解读
EM算法 - 期望极大算法