1. 定义
最大期望算法(Expectation-maximizationalgorithm,又译为期望最大化算法),是在概率模型中寻找参数最大似然估计或者最大后验估计的算法,其中概率模型依赖于无法观测的隐性变量。2. 算法步骤
最大期望算法经过两个步骤交替进行计算:
第一步是计算期望(E),利用对隐藏变量的现有估计值,计算其最大似然估计值;
第二步是最大化(M),最大化在E步上求得的最大似然值来计算参数的值。M步上找到的参数估计值被用于下一个E步计算中,这个过程不断交替进行。3. 什么是似然函数?
在数理统计学中,似然函数是一种关于统计模型中的参数的函数,表示模型参数中的似然性。“似然性”与“或然性”或“概率”意思相近,都是指某种事件发生的可能性。
而极大似然就相当于最大可能的意思。
比如你一位同学和一位猎人一起外出打猎,一只野兔从前方窜过。只听一声枪响,野兔应声到下,如果要你推测,这一发命中的子弹是谁打的?你就会想,只发一枪便打中,由于猎人命中的概率一般大于你那位同学命中的概率,从而推断出这一枪应该是猎人射中的。
这个例子所作的推断就体现了最大似然法的基本思想。
多数情况下我们是根据已知条件来推算结果,而最大似然估计是已经知道了结果,然后寻求使该结果出现的可能性最大的条件,以此作为估计值。
概率是根据条件推测结果,而似然则是根据结果反推条件。在这种意义上,似然函数可以理解为条件概率的逆反。
假定已知某个参数B时,推测事件A会发生的概率写作:
通过贝叶斯公式,可以得出
求极大似然函数估计值的一般步骤:
[1]写出似然函数;
[2]对似然函数取对数,并整理;
[3]求导数,令导数为0,得到似然方程;
[4]解似然方程,得到的参数即为所求。
4. EM算法的经典案例:抛硬币
两枚硬币A和B,假定随机抛掷后正面朝上概率分别为PA,PB。为了估计这两个硬币朝上的概率,轮流抛硬币A和B,每一轮都连续抛5次,总共5轮:
硬币A被抛了15次,在第一轮、第三轮、第五轮分别出现了3次正、1次正、2次正,所以很容易估计出PA,类似的,PB也很容易计算出来,如下:
PA =(3+1+2)/ 15 = 0.4
PB=(2+3)/10 = 0.5
问题来了,如果我们不知道抛的硬币是A还是B呢(即硬币种类是隐变量),然后再轮流抛五轮,得到如下结果:
问题变得有意思了。现在我们的目标没变,还是估计PA和PB,需要怎么做呢?
显然,此时我们多了一个硬币种类的隐变量,设为z,可以把它认为是一个5维的向量(z1,z2,z3,z4,z5),代表每次投掷时所使用的硬币,比如z1,就代表第一轮投掷时使用的硬币是A还是B。
但是,这个变量z不知道,就无法去估计PA和PB,所以,我们必须先估计出z,然后才能进一步估计PA和PB。可要估计z,我们又得知道PA和PB,这样我们才能用极大似然概率法则去估计z。
答案就是先随机初始化一个PA和PB,用它来估计z,然后基于z,还是按照最大似然概率法则去估计新的PA和PB,不断迭代,直至收敛。
我们不妨这样,先随便给PA和PB赋一个值,比如:
硬币A正面朝上的概率PA = 0.2
硬币B正面朝上的概率PB = 0.7
然后,我们看看第一轮抛掷最可能是哪个硬币。
如果是硬币A,得出3正2反的概率为0.2*0.2*0.2*0.8*0.8 = 0.00512
如果是硬币B,得出3正2反的概率为0.7*0.7*0.7*0.3*0.3=0.03087
然后依次求出其他4轮中的相应概率。做成表格如下(标粗表示其概率更大):
按照最大似然法则:
第1轮中最有可能的是硬币B,第2轮中最有可能的是硬币A,第3轮中最有可能的是硬币A,第4轮中最有可能的是硬币B,第5轮中最有可能的是硬币A。
我们就把概率更大,即更可能是A的,即第2轮、第3轮、第5轮出现正的次数2、1、2相加,除以A被抛的总次数15(A抛了三轮,每轮5次),作为z的估计值,B的计算方法类似。然后我们便可以按照最大似然概率法则来估计新的PA和PB。
PA =(2+1+2)/15 = 0.33
PB =(3+3)/10 = 0.6
设想我们是全知的神,知道每轮抛掷时的硬币就是如本文本节开头标示的那样,那么,PA和PB的最大似然估计就是0.4和0.5(下文中将这两个值称为PA和PB的真实值)。那么对比下我们初始化的PA和PB和新估计出的PA和PB。
就这样,不断迭代不断接近真实值,这就是EM算法的奇妙之处。
可以期待,我们继续按照上面的思路,用估计出的PA和PB再来估计z,再用z来估计新的PA和PB,反复迭代下去,找到保持不变的值。我们就找到了PA和PB的最大似然估计。
假设我们有一个样本集{x(1),…,x(m)},包含m个独立的样本,现在我们想找到每个样本隐含的类别z,能使得p(x,z)最大。p(x,z)的极大似然估计如下:
第一步是对极大似然取对数,第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求一般比较困难,因为有隐藏变量z存在,但是一般确定了z后,求解就容易了。
对于参数估计,我们本质上还是想获得一个使似然函数最大化的那个参数θ,现在与极大似然不同的只是似然函数式中多了一个未知的变量z,见下式(1)。也就是说我们的目标是找到适合的θ和z,以让L(θ)最大。那我们也许会想,你就是多了一个未知的变量而已啊,我也可以分别对未知的θ和z分别求偏导,再令其等于0,求解出来不也一样吗?
本质上我们是需要最大化(1)式,也就是似然函数,我们把分子分母同乘以一个相等的函数(即隐变量Z的概率分布Qi(z(i)),其概率之和等于1,即,即得到上图中的(2)式,我们通过Jensen不等式可得到(3)式,此时(3)式变成了“对数的和”,如此求导就容易了。
当我们在求(2)式的极大值时,(2)式不太好计算,我们便想(2)式大于某个方便计算的下界(3)式,而我们尽可能让下界(3)式最大化,直到(3)式的计算结果等于(2)式,便也就间接求得了(2)式的极大值;
关于Jensen不等式,定义如下:如果f是凸函数,X是随机变量,那么:E[f(X)]>=f(E[X]),通俗的说法是函数的期望大于等于期望的函数。特别地,如果f是严格凸函数,当且仅当P(X = EX) = 1,即X是常量时,上式取等号,即E[f(X)] = f(EX)。
5. EM算法应用
关于EM算法的应用其实很多,最广泛的就是GMM混合高斯模型、聚类、HMM等等。如下是一个在Python 3环境使用EM算法求解GMM的编程实现。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
#构建测试数据
N = 200; pi1 = np.array([0.6, 0.3, 0.1])
mu1 = np.array([[0,4], [-2,0], [3,-3]])
sigma1 = np.array([[[3,0],[0,0.5]], [[1,0],[0,2]], [[.5,0],[0,.5]]])
gen = [np.random.multivariate_normal(mu, sigma, int(pi*N)) for mu,sigma, pi in zip(mu1, sigma1, pi1)]
X = np.concatenate(gen)
#初始化: mu, sigma, pi =均值,协方差矩阵,混合系数
theta = {}; param = {}
theta['pi'] = [1/3, 1/3, 1/3] #均匀初始化
theta['mu'] = np.random.random((3, 2)) #随机初始化
theta['sigma'] = np.array([np.eye(2)]*3) #初始化为单位正定矩阵
param['k'] = len(pi1); param['N'] = X.shape[0]; param['dim'] =X.shape[1]
#定义函数
def GMM_component(X, theta, c):
'''
由联合正态分布计算GMM的单个成员
'''
returntheta['pi'][c]*multivariate_normal(theta['mu'][c], theta['sigma'][c,...]).pdf(X)
def E_step(theta, param):
'''
E步:更新隐变量概率分布q(Z)。
'''
q = np.zeros((param['k'],param['N']))
for i in range(param['k']):
q[i, :] = GMM_component(X,theta, i)
q /= q.sum(axis=0)
return q
def M_step(X, q, theta, param):
'''
M步:使用q(Z)更新GMM参数。
'''
pi_temp = q.sum(axis=1);pi_temp /= param['N'] #计算pi
mu_temp = q.dot(X); mu_temp /=q.sum(axis=1)[:, None] #计算mu
sigma_temp =np.zeros((param['k'], param['dim'], param['dim']))
for i in range(param['k']):
ys = X - mu_temp[i, :]
sigma_temp[i] =np.sum(q[i, :, None, None]*np.matmul(ys[..., None], ys[:, None, :]), axis=0)
sigma_temp /= np.sum(q,axis=1)[:, None, None] #计算sigma
theta['pi'] = pi_temp; theta['mu']= mu_temp; theta['sigma'] = sigma_temp
return theta
def likelihood(X, theta):
'''
计算GMM的对数似然。
'''
ll = 0
for i in range(param['k']):
ll += GMM_component(X,theta, i)
ll = np.log(ll).sum()
return ll
def EM_GMM(X, theta, param, eps=1e-5, max_iter=1000):
'''
高斯混合模型的EM算法求解
theta: GMM模型参数; param:其它系数; eps:计算精度; max_iter:最大迭代次数
返回对数似然和参数theta,theta是包含pi、mu、sigma的Python字典
'''
for i in range(max_iter):
ll_old = 0
# E-step
q = E_step(theta, param)
# M-step
theta = M_step(X, q,theta, param)
ll_new = likelihood(X,theta)
if np.abs(ll_new - ll_old)< eps:
break;
else:
ll_old = ll_new
return ll_new, theta
# EM算法求解GMM,最大迭代次数为1e5
ll, theta2 = EM_GMM(X, theta, param, max_iter=10000)
#由theta计算联合正态分布的概率密度
L = 100; Xlim = [-6, 6]; Ylim = [-6, 6]
XGrid, YGrid = np.meshgrid(np.linspace(Xlim[0], Xlim[1], L),np.linspace(Ylim[0], Ylim[1], L))
Xout = np.vstack([XGrid.ravel(), YGrid.ravel()]).T
MVN = np.zeros(L*L)
for i in range(param['k']):
MVN += GMM_component(Xout,theta, i)
MVN = MVN.reshape((L, L))
#绘制结果
plt.plot(X[:, 0], X[:, 1], 'x', c='gray', zorder=1)
plt.contour(XGrid, YGrid, MVN, 5, colors=('k',), linewidths=(2,))
参考文献:July在其CSDN博客上的EM笔记《如何通俗理解EM算法》。