複製鏈接
請複製以下鏈接發送給好友
https://baike.baidu.hk/item/最大期望算法/10180861
複製
複製成功

最大期望算法

鎖定
最大期望算法（Expectation-Maximization algorithm, EM），或Dempster-Laird-Rubin算法 [1]  ，是一類通過迭代進行極大似然估計（Maximum Likelihood Estimation, MLE）的優化算法 [2]  ，通常作為牛頓迭代法（Newton-Raphson method）的替代用於對包含隱變量（latent variable）或缺失數據（incomplete-data）的概率模型進行參數估計 [2-3]  。
EM算法的標準計算框架由E步（Expectation-step）和M步（Maximization step）交替組成，算法的收斂性可以確保迭代至少逼近局部極大值 [4]  。EM算法是MM算法（Minorize-Maximization algorithm）的特例之一，有多個改進版本，包括使用了貝葉斯推斷的EM算法、EM梯度算法、廣義EM算法等 [2]  。
由於迭代規則容易實現並可以靈活考慮隱變量 [3]  ，EM算法被廣泛應用於處理數據的缺測值 [1-2]  ，以及很多機器學習（machine learning）算法，包括高斯混合模型（Gaussian Mixture Model, GMM） [5]  和隱馬爾可夫模型（Hidden Markov Model, HMM） [6]  的參數估計。
中文名
最大期望算法
外文名
Expectation Maximization algorithm, EM
類    型
優化算法
提出者
Arthur Dempster，Nan Laird，
自定義
Donald Rubin 等
提出時間
1977年
應    用
數據分析，機器學習

目錄

  1. 1 歷史
  2. 2 理論
  3. 3 算法
  1. ▪ 標準算法
  2. ▪ 改進算法
  3. 4 性質
  1. 5 應用

最大期望算法歷史

對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）理論的優化算法。給定相互獨立的觀測數據
，和包含隱變量
、參數
的概率模型
，根據MLE理論，
的最優單點估計在模型的似然取極大值時給出：
。考慮隱變量，模型的似然有如下展開 [15]  [3]  ：
隱變量可以表示缺失數據，或概率模型中任何無法直接觀測的隨機變量，式中第一行是隱變量為連續變量的情形，第二行是隱變量為離散變量的情形，積分/求和的部分也被稱為
的聯合似然（joint liklihood）。不失一般性，這裏按離散變量為例進行説明。由MLE的一般方法，對上式取自然對數後可得 [3]  ：
上述展開考慮了觀測數據的相互獨立性。引入與隱變量有關的概率分佈
，即隱分佈（可認為隱分佈是隱變量對觀測數據的後驗，參見標準算法的E步推導），由Jensen不等式，觀測數據的對數似然有如下不等關係 [3]  ：
當
使不等式右側取全局極大值時，所得到的
至少使不等式左側取局部極大值。因此，將不等式右側表示為
後，EM算法有如下求解目標 [3]  ：
式中的
等效於MM算法（Minorize-Maximization algorithm）中的代理函數（surrogate function），是MLE優化問題的下限，EM算法通過最大化代理函數逼近對數似然的極大值。

最大期望算法算法

最大期望算法標準算法

計算框架
EM的計算框架：對數似然（藍），E步（綠），M步（實心點） EM的計算框架：對數似然（藍），E步（綠），M步（實心點） [16]
EM標準算法是一組迭代計算，迭代分為兩部分，即E步和M步，其中E步“固定”前一次迭代的
，求解
使
取極大值；M步使用
求解
使
取極大值。EM算法在初始化模型參數後開始迭代，迭代中E步和M步交替進行。由於EM算法的收斂性僅能確保局部最優，而不是全局最優 [3-4]  。因此通常對EM算法進行隨機初始化並多次運行，選擇對數似然最大的迭代輸出結果 [3]  。以下給出EM算法E步和M步的推導。
1. E步（Expectation-step, E-step）
由EM算法的求解目標可知，E步有如下優化問題 [3]  ：
考慮先前的不等關係，這裏首先對
進行展開 [3]  ：
注意到，推導上式時考慮了：
。由貝葉斯定理（Bayes' theorem），上式可化為 [3]  ：
E步：優化代理損失（左），原優化目標（右） E步：優化代理損失（左），原優化目標（右） [15]
式中
為Kullback-Leibler散度（Kullback-Leibler divergence, KL）或相對熵（relative entropy），
表示吉布斯自由能（Gibbs free energy），即由Jensen不等式得到的代理函數等價於隱分佈的自由能。求解
的極大值等價於求解隱分佈自由能的極大值，即隱分佈對隱變量後驗
的KL散度的極小值。由KL散度的性質可知，其極小值在兩個概率分佈相等時取得，因此當
時，
取極大值，對EM算法的第
次迭代，E步有如下計算 [3]  [17]  ：
2. M步（Maximization step, M-step）
在E步的基礎上，M步求解模型參數使
取極大值。該極值問題的必要條件是
[3]  ：
式中
表示聯合似然
對隱分佈
的數學期望。在
為凸函數時（例如隱變量和觀測服從指數族分佈），上述推導也是充分的 [3]  。由此得到M步的計算：
計算步驟
將統計模型帶入EM算法的計算框架即可得到其計算步驟。這裏以高斯混合模型（Gaussian Mixture Model, GMM）為例進行説明。由GMM的一般定義可知，其似然和參數有如下表示 [3]  [5]  ：
根據學習數據的維度，式中
表示均值為
，方差/協方差為
的正態分佈/聯合正態分佈。
為正態分佈的混合比例，
為參與混合的分佈總數。定義與觀測數據有關的隱變量：
，令隱分佈
表示GMM聚類的軟指定（soft assignment），即每個數據來源於第
個分佈的概率，則隱變量有離散取值
。
將上述內容帶入EM算法的計算框架後，E步有如下展開 [17]  ：
GMM中有：
，因此E步的計算步驟為：
M步通過E步輸出的隱變量後驗計算模型參數，在GMM中，M步計算框架的優化問題有如下表示 [17]  ：
不失一般性，帶入單變量正態分佈的解析形式後對模型參數求偏導數可得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算法僅能給出模型參數
的單點估計，引入貝葉斯推斷方法後，EM算法能夠給出模型參數的後驗（posterior）分佈避免過度擬合，其中常見的例子是極大後驗估計（Maximum A Posteriori estimation, MAP）的EM算法 [2]  [18]  。MAP-EM在標準EM算法的基礎上引入了模型參數的先驗（prior）：
，此時MAP-EM的優化目標由模型的似然轉變為後驗，其離散形式可表示為 [18]  ：
類比標準EM算法，考慮隱分佈
後，由Jensen不等式可得到對數後驗的代理函數，即隱變量的自由能 [18]  ：
由此可得MAP-EM的迭代步驟：
MAP-EM在Dempster et al. (1977)就已被提出 [1]  ，但不同於標準EM，MAP-EM的隱分佈
是隱變量和模型參數的聯合分佈，其對應的隱變量後驗
往往沒有解析形式。在貝葉斯體系下，求解該隱變量後驗的方法包括馬爾可夫鏈蒙特卡羅（Markov Chain Monte Carlo, MCMC）和變分貝葉斯估計（Variational Bayesian Inference, VBI），對前者，可證明由MCMC求解的MAP-EM等價於吉布斯採樣（Gibbs sampling）算法 [19]  ；對後者，由VBI求解的MAP-EM被稱為變分貝葉斯EM算法（Variational Bayesian EM algorithm, VBEM） [18]  [3]  。
VBEM使用平均場理論（Mean Field Theory, MFT）將隱分佈近似為其在各個維度上分佈的乘積：
並由此得到以下迭代步驟 [2]  ：
使用VBEM的常見例子是語言建模問題中的隱含狄利克雷分佈（latent dirichlet allocation） [3]  。
EM梯度算法（EM gradient algorithm）和廣義EM算法（Generalized EM algorithm, GEM）
參見：廣義期望最大算法
EM算法的M步通過計算偏導數求解代理函數的極大值，EM梯度算法（EM gradient algorithm）將該過程替換為牛頓迭代法（Newton-Raphson method）以加速迭代收斂 [2]  [20]  。更進一步地，當代理函數
不是凸函數或無法有效地對
求解極大值時，可以使用廣義EM算法（GEM）。GEM有兩種實現方式，一是在M步使用非線性優化策略，例如共軛梯度算法（conjugate gradients algorithm） [21]  ，二是將原M步的求導計算分解為多個條件極值問題逐個計算模型參數，後者也被稱為最大條件期望算法（Expectation Conditional Maximization algorithm, ECM） [15]  。
EM算法的E步也可以按ECM的方法分解為條件極值問題，由先前推導可知，E步的優化問題僅有一個全局極大值，即隱分佈
，因此在E步將MLE的優化目標：聯合似然
對觀測樣本按因子展開並對每個展開分別使用EM算法，可以得到同樣的優化結果。對於M步，如果觀測樣本來自指數族分佈，則M步也可以在每次迭代僅對有限個樣本的展開進行。在指數族問題中使用EM算法的上述推廣，可以避免在迭代中反覆處理整個觀測樣本集，降低計算開銷 [15]  [22]  。
α-EM算法（α-EM algorithm）
α-EM算法是對標準算法的隱變量概率分佈引入權重係數
的改進版本。標準的EM算法是α-EM算法在
時的特例。給定恰當的超參數
，α-EM能夠比標準EM算法更快收斂。有研究將α-EM算法應用於神經網絡的概率學習和隱馬爾可夫模型的參數估計 [23-24]  。

最大期望算法性質

收斂性與穩定性：EM算法必然收斂於對數似然的局部極大值或鞍點（saddle point），其證明考慮如下不等關係 [3]  ：
由上式可知EM算法得到的對數似然是單調遞增的，即從
次迭代到
次迭代，EM算法至少能維持當前的優化結果，不會向極大值的相反方向運動，因此EM算法具有數值穩定性（numerical stablility）。上述不等關係也被用於EM算法迭代終止的判定，給定計算精度
，當
時迭代結束。EM算法收斂性的具體證明參見Wu (1983) [4]  。
計算複雜度：在E步具有解析形式時，EM算法是一個計算複雜度和存儲開銷都很低的算法，可以在很小的計算資源下完成計算 [2]  。在E步不具有解析形式或使用MAP-EM時，EM算法需要結合其它數值方法，例如變分貝葉斯估計或MCMC對隱變量的後驗分佈進行估計，此時的計算開銷取決於問題本身 [2]  [3]  。
與其它算法的比較：相比於梯度算法，例如牛頓迭代法和隨機梯度下降（Stochastic Gradient Descent, SGD），EM算法的優勢是其求解框架可以加入求解目標的額外約束，例如在高斯混合模型的例子中，EM算法在求解協方差時可以確保每次迭代的結果都是正定矩陣 [3]  。EM算法的不足在於其會陷入局部最優，在高維數據的問題中，局部最優和全局最優可能有很大差異 [2]  。

最大期望算法應用

EM算法及其改進版本被用於機器學習算法的參數求解，常見的例子包括高斯混合模型（Gaussian Mixture Model, GMM） [5]  、概率主成份分析（probabilistic Principal Component Analysis） [25]  、隱馬爾可夫模型（Hidden Markov Model, HMM） [6]  等非監督學習算法。EM算法可以給出隱變量，即缺失數據的後驗
，因此在缺失數據問題（incomplete-data probelm）中有應用 [1]  [2]  。
參考資料
  • 1.    Dempster, A.P., Laird, N.M. and Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. Journal of the royal statistical society. Series B (methodological), pp.1-38.
  • 2.    McLachlan, G. and Krishnan, T.．The EM algorithm and extensions (Vol. 382)：John Wiley & Sons.，2007：pp.1-73, 276-277
  • 3.    Polykovskiy, D. and Novikov, A., Bayesian Methods for Machine Learning  ．Coursera and National Research University Higher School of Economics．2017[引用日期2018-12-12]
  • 4.    Wu, C.J., 1983. On the convergence properties of the EM algorithm. The Annals of statistics, pp.95-103.
  • 5.    Xu, L. and Jordan, M.I., 1996. On convergence properties of the EM algorithm for Gaussian mixtures. Neural computation, 8(1), pp.129-151.
  • 6.    Bilmes, J.A., 1998. A gentle tutorial of the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov models. International Computer Science Institute, 4(510), p.126.
  • 7.    Newcomb, S., 1886. A generalized theory of the combination of observations so as to obtain the best result. American journal of Mathematics, pp.343-366.
  • 8.    Dietz, K., 1997. Introduction to McKendrick (1926) Applications of mathematics to medical problems. In Breakthroughs in Statistics (pp. 17-57). Springer, New York, NY.
  • 9.    Healy, M. and Westmacott, M., 1956. Missing values in experiments analysed on automatic computers. Applied statistics, pp.203-206.
  • 10.    Blight, B.J.N., 1970. Estimation from a censored sample for the exponential family. Biometrika, 57(2), pp.389-395.
  • 11.    Sundberg, R., 1972. Maximum Likelihood Theory and Applications for Distributions Generated when Observing a Function an Exponential Family Variable. Institute of Mathematical Statics, Stockholm University.
  • 12.    Sundberg, R., 1974. Maximum likelihood theory for incomplete data from an exponential family. Scandinavian Journal of Statistics, pp.49-58.
  • 13.    Baum, L.E. and Petrie, T., 1966. Statistical inference for probabilistic functions of finite state Markov chains. The annals of mathematical statistics, 37(6), pp.1554-1563.
  • 14.    Baum, L.E. and Eagon, J.A., 1967. An inequality with applications to statistical estimation for probabilistic functions of Markov processes and to a model for ecology. Bulletin of the American Mathematical Society, 73(3), pp.360-363.
  • 15.    Bishop, C.M.．Pattern recognition and machine learning. In information science and statistics．New York, NJ：Springer，2006：450-455
  • 16.    Do, C.B. and Batzoglou, S., 2008. What is the expectation maximization algorithm?. Nature biotechnology, 26(8), supplementary note, p.897.
  • 17.    Expecation Maximization, STA-663: Computational Statistics and Statistical Computing  ．Duke University．2016[引用日期2018-12-13]
  • 18.    Beal, M.J., 2003. Variational algorithms for approximate Bayesian inference. Doctor of Philosophy Thesis. University of London.
  • 19.    Mitchell, T., Gibbs and Metropolis sampling (MCMC methods) and relations of Gibbs to EM. In: 10-701 Machine Learning  ．School of Computer Science, Carnegie Mellon University．2011[引用日期2019-01-22]
  • 20.    Lindstrom, M.J. and Bates, D.M., 1988. Newton—Raphson and EM algorithms for linear mixed-effects models for repeated-measures data. Journal of the American Statistical Association, 83(404), pp.1014-1022.
  • 21.    Meng, X.L. and Rubin, D.B., 1993. Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika, 80(2), pp.267-278.
  • 22.    Neal, R.M. and Hinton, G.E., 1998. A view of the EM algorithm that justifies incremental, sparse, and other variants. In Learning in graphical models (pp. 355-368). Springer, Dordrecht.
  • 23.    Matsuyama, Y., 1997. The α-EM algorithm: A block connectable generalized leaning tool for neural networks. In International Work-Conference on Artificial Neural Networks (pp. 483-492). Springer, Berlin, Heidelberg.
  • 24.    Matsuyama, Y., 2011. Hidden Markov model estimation based on alpha-EM algorithm: Discrete and continuous alpha-HMMs. In Neural Networks (IJCNN), The 2011 International Joint Conference on (pp. 808-816). IEEE.
  • 25.    Tipping, M.E. and Bishop, C.M., 1999. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3), pp.611-622.
展開全部 收起