机器学习之EM算法
- 1、EM算法介绍
- 2、EM算法数学原理
- 3、算法及Python实现
- 4、小结
1、EM算法介绍
EM算法引入:
首先使用经典的三硬币模型引入EM算法。
三硬币模型
有A,B,C三枚硬币,单次投掷出现正面的概率分别为π、p、q。利用这三枚硬币进行如下实验:
1、先投掷A,若出现正面则投掷B,否则投掷C
2、记录第二次投掷的硬币出现的结果,正面记作1,反面记作0
独立重复步骤1和步骤2十次,观察结果如下:
1 1 0 1 0 0 1 0 1 1
假设只能观测到投掷硬币结果,不能观测投掷硬币的过程,求π、p、q,即三硬币模型的参数。
记模型参数θ=(π,p,q),无法观测的第二次投掷的硬币为隐变量z,y为可以观测到第二次投掷的结果0或1,则观测数据的释然函数为
这个式子表示两个事件,第一个事件选取那枚看不到的硬币,第二个事件利用这枚硬币进行一次投掷。利用硬币的结果只可能是0或1这个特性,可以将式子展开为:
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/4ba9841f4899fb8df62f0bd33e31dfea.jpeg)
y的观测序列给定了,那么就是找到一个模型参数使得这个序列的概率(似然函数的值)最大,也就是求模型参数的极大似然估计:
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/83daa07a914b2aae73c5ccd62bc953d6.jpeg)
2、EM算法的数学原理
EM算法是求解三硬币这类问题的一种迭代算法,它有3步:
初始化:选取模型参数的初值:,循环如下两步迭代
E步:计算在单签迭代的模型参数下,观测数据y来自硬币B的概率:
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/953da16b4f73b10e6e75f10307e7cf87.jpeg)
这个式子中分子代表选定B并进行一次投掷实验,分母代表选定B或C并进行一次投掷实验,两个一除就得到实验结果来自B的概率。
M步:估算下一个迭代的新的模型估算值:
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/073cbadd01ea7962adb30db43ca897c9.jpeg)
整体思路是把n个{实验结果来自B的概率}求和得到期望,平均后,得到B出现正面的似然估计,同理有p和q。
重复迭代,直到收敛为止。
这个模型中,观测数据Y和隐数据Z组合在一起称为完全数据,单独的观测数据Y称为不完全数据。在隐数据未知的情况,无法直接估计Y的概率分布。但当模型概率给定时,就可以估计Y的条件概率分布了。
3、算法及Python实现
EM算法的标准定义
输入:观测变量数据Y,隐变量数据Z,联合概率分布P(Y,Z|θ),条件分布P(Z|Y,θ);
输出:模型参数θ
(1)选择参数的初始值
θ(0)
θ
(
0
)
,开始迭代;
(2)E步:记
θ(i)
θ
(
i
)
为第i次迭代参数θ的估计值,在第i+1次迭代的E步,计算
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/c4d41077e9727e79e8ed2ceb64571382.jpeg)
这里,P(Z|Y, θ(i) θ ( i ) )是在给定观测数据Y和当前参数估计 θ(i) θ ( i ) 下隐变量数据Z的条件概率分布;
(3)M步:求使Q(θ, θ(i) θ ( i ) )极大化的θ,确定第i+1次迭代的参数的估计值 θ(i) θ ( i )
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/466357575253b433cdf59f923f667922.jpeg)
(4)重复第(2)和(3)步,直到收敛。
式子
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/f7f118d3fa1a2017020cd1709c4210de.jpeg)
定义(Q函数)完全数据的对数似然函数 logP(Y,Z|θ) l o g P ( Y , Z | θ ) 关于在给定观测数据Y和当前参数 θ(i) θ ( i ) 下对未观测数据Z的条件概率分布P(Z|Y, θ(i) θ ( i ) )的期望称为Q函数,即
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/11da05f2fbc1d8348447a8b8918e80b9.jpeg)
下面关于EM算法作几点说明:
步骤(1)参数的初值可以任意选择,但需注意EM算法对初值是敏感的。
步骤(2)E步求Q(θ, θ(i) θ ( i ) )。Q函数式中Z是未观测数据,Y是观测数据。注意,Q(θ, θ(i) θ ( i ) )的第1个变元表示要极大化的参数,第2个变元表示参数的当前估计值。每次迭代实际在求Q函数及其极大。
步骤(3)M步求Q(θ, θ(i) θ ( i ) )的极大化,得到 θ(i+1) θ ( i + 1 ) ,完成一次迭代 θ(i)→θ(i+1) θ ( i ) → θ ( i + 1 ) 。后面将证明每次迭代使似然函数增大或达到局部极值。
步骤(4)给出停止迭代的条件,一般是对较小的正数
ε1,ε2
ε
1
,
ε
2
,若满足
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/7df5318edf1158725774fbf0f6262abe.jpeg)
EM算法的简单应用
假设有两枚硬币A、B,以相同的概率随机选择一个硬币,进行如下的抛硬币实验:共做5次实验,每次实验独立的抛十次,结果如图中a所示,例如某次实验产生了H、T、T、T、H、H、T、H、T、H,H代表正面朝上。
但是由于在记录实验时没能记录每次实验选择的是硬币A还是B,我们无法观测实验数据中选择的硬币是哪个,问在这种情况下两个硬币正面出现的概率?
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/a9723e47561c3fc79852d0cf2f0b039d.jpeg)
下面采用EM算法解决该问题
Python代码如下
import scipy.stats as stats
import math
import numpy as np
# 硬币投掷结果观测序列
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]])
def em_single(priors,observations):
"""
EM算法单次迭代过程
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}} #AB硬币统计正反面次数,H正面,T反面
theta_A = priors[0]
theta_B = priors[1]
# E Step
for observation in observations:
len_observation = len(observation)
num_heads = observation.sum()
num_tails = len_observation - num_heads
contribution_A = stats.binom.pmf(num_heads,len_observation,theta_A)
contribution_B = stats.binom.pmf(num_heads,len_observation,theta_B)
weight_A = contribution_A / (contribution_A + contribution_B)
weight_B = contribution_B / (contribution_A + contribution_B)
#更新在当前参数下A、B硬币的正反面次数
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 step
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]
def em(observation,prior,tol=1e-6,iterations=10000):
"""
EM算法
:param observation: 观测数据
:param prior: 模型初值
:param tol: 迭代结束阈值
:param iterations: 最大迭代次数
:return: 局部最优的模型参数
"""
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]
[prob_A,prob_B],iteration = em(observations,[0.6,0.5])
print("The probability of A is %f,the probability of B is %f"%(prob_A,prob_B))
print('The iteration is %d'%iteration)
输出结果如下:
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/db9c24a5da1a23ca16f8d96af4591830.png)
4、小结
EM算法寻找模型参数的目标(或称标准)是找到的参数使观测数据的似然函数最大,一般用对数似然函数取代似然函数,这样可以把连乘变为累加,方便优化,也就是极大化
![这里写图片描述](https://i-blog.csdnimg.cn/blog_migrate/30a783586020522568e7fa83cd51bbc0.jpeg)
这个式子里面有未知的隐变量Z,无法直接优化。
但是如同在“EM算法简单理解”中看到那样,给定模型参数,就可以估计Y的条件概率(后验概率,已经有Z这个结果,求原因Y的概率)。所以我们就挑一个模型参数的初值,也就是EM算法的第1步。
有了初值,就可以代入似然函数得到一个值,但这个值不一定是最大的,我们想要更大,所以需要调整参数,这也是EM算法为什么要迭代的原因。
事实上,EM算法是通过迭代逐步近似极大化似然函数的。假设在第i次迭代后θ的估计值是
θ(i)
θ
(
i
)
。我们希望新估计值θ能使L(θ)增加,即L(θ)>L(
θ(i)
θ
(
i
)
),并逐步达到极大值。为此,考虑两者之差:
利用Jensen不等式(Jensen inequality)
得到其下界:
第一行的
P(Y|Z,θ(i))
P
(
Y
|
Z
,
θ
(
i
)
)
是人为加上去的,先乘以这一项再除以这一项结果仍然是L(θ),然后第二行就利用了Jensen不等式,将log运算符移入求和项中。但
P(Y|Z,θ(i))
P
(
Y
|
Z
,
θ
(
i
)
)
后面变成了
P(Z|Y,θ(i))
P
(
Z
|
Y
,
θ
(
i
)
)
,这里可能是李航博士的笔误导致的,第一行时分子分母应该同时乘以
P(Z|Y,θ(i))
P
(
Z
|
Y
,
θ
(
i
)
)
的。参考普林斯顿大学的讲义《COS 424- Interacting with Data.pdf》