最大期望算法
栏目：行业资讯 发布时间：2024-08-26
播报编辑对EM算法的研究起源于统计学的误差分析（erroranalysis）问题。1886年，美国数学家SimonNewcomb在使用高斯混合模型（GaussianMixtureModel,GMM）解释观测误差的长尾效应时提出了类似EM算法的迭代求解技术[7]。在极大似然估计（MaximumLik
对EM算法的研究起源于统计学的误差分析（error analysis）问题。1886年，美国数学家Simon Newcomb在使用高斯混合模型（Gaussian Mixture Model, GMM）解释观测误差的长尾效应时提出了类似EM算法的迭代求解技术 [7]。在极大似然估计（Maximum Likelihood Estimation, MLE）方法出现后，英国学者Anderson McKendrick在1926年发展了Newcomb的理论并在医学样本中进行了应用 [8]。1956年，Michael Healy和Michael Westmacott提出了统计学试验中估计缺失数据的迭代方法 [9]，该方法被认为是EM算法的一个特例 [2]。1970年，B. J. N. Blight使用MLE对指数族分布的I型删失数据（Type I censored data）进行了讨论 [10]。Rolf Sundberg在1971至1974年进一步发展了指数族分布样本的MLE并给出了迭代计算的完整推导 [11-12]。
EM算法的正式提出来自美国数学家Arthur Dempster、Nan Laird和Donald Rubin，其在1977年发表的研究对先前出现的作为特例的EM算法进行了总结并给出了标准算法的计算步骤，EM算法也由此被称为Dempster-Laird-Rubin算法 [1-2]。1983年，美国数学家吴建福（C.F. Jeff Wu）给出了EM算法在指数族分布以外的收敛性证明 [4]。
此外，在二十世纪60-70年代对隐马尔可夫模型（Hidden Markov Model, HMM）的研究中，Leonard E. Baum提出的基于MLE的HMM参数估计方法，即Baum-Welch算法（Baum-Welch algorithm）也是EM算法的特例之一 [6] [13-14]。
EM算法是基于极大似然估计（Maximum Likelihood Estimation, MLE）理论的优化算法。给定相互独立的观测数据
计算框架
EM的计算框架：对数似然（蓝），E步（绿），M步（实心点） [16]
1. E步（Expectation-step, E-step）
由EM算法的求解目标可知，E步有如下优化问题 [3]：
E步：优化代理损失（左），原优化目标（右） [15]
2. M步（Maximization step, M-step）
在E步的基础上，M步求解模型参数使
计算步骤
将统计模型带入EM算法的计算框架即可得到其计算步骤。这里以高斯混合模型（Gaussian Mixture Model, GMM）为例进行说明。由GMM的一般定义可知，其似然和参数有如下表示 [3] [5]：
将上述内容带入EM算法的计算框架后，E步有如下展开 [17]：
M步通过E步输出的隐变量后验计算模型参数，在GMM中，M步计算框架的优化问题有如下表示 [17]：
根据以上计算步骤，这里给出一个在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的单个成员
'''
return?theta['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字典
'''
ll_old?=?0
for?i?in?range(max_iter):
#?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,))
基于贝叶斯推断（Bayesian inference）的EM算法
参见：变分贝叶斯EM
在MLE理论下，EM算法仅能给出模型参数
MAP-EM在Dempster et al. (1977)就已被提出 [1]，但不同于标准EM，MAP-EM的隐分布
EM梯度算法（EM gradient algorithm）和广义EM算法（Generalized EM algorithm, GEM）
参见：广义期望最大算法
EM算法的M步通过计算偏导数求解代理函数的极大值，EM梯度算法（EM gradient algorithm）将该过程替换为牛顿迭代法（Newton-Raphson method）以加速迭代收敛 [2] [20]。更进一步地，当代理函数
EM算法的E步也可以按ECM的方法分解为条件极值问题，由先前推导可知，E步的优化问题仅有一个全局极大值，即隐分布
α-EM算法（α-EM algorithm）
α-EM算法是对标准算法的隐变量概率分布引入权重系数
计算复杂度：在E步具有解析形式时，EM算法是一个计算复杂度和存储开销都很低的算法，可以在很小的计算资源下完成计算 [2]。在E步不具有解析形式或使用MAP-EM时，EM算法需要结合其它数值方法，例如变分贝叶斯估计或MCMC对隐变量的后验分布进行估计，此时的计算开销取决于问题本身 [2-3]。
与其它算法的比较：相比于梯度算法，例如牛顿迭代法和随机梯度下降（Stochastic Gradient Descent, SGD），EM算法的优势是其求解框架可以加入求解目标的额外约束，例如在高斯混合模型的例子中，EM算法在求解协方差时可以确保每次迭代的结果都是正定矩阵 [3]。EM算法的不足在于其会陷入局部最优，在高维数据的问题中，局部最优和全局最优可能有很大差异 [2]。
