EM算法也稱期望最大化（Expectation-Maximum, EM）算法，它是一個基礎算法，是很多機器學習領域算法的基礎，比如隱式馬爾科夫算法（HMM）， LDA主題模型的變分推斷等。本文就對EM算法的原理做一個總結。

?

我們經常會從樣本觀察數據中，找出樣本的模型參數。 最常用的方法就是極大化模型分布的對數似然函數。但是在一些情況下，我們得到的觀察數據有未觀察到的隱含數據，此時我們未知的有隱含數據和模型參數，因而無法直接用極大化對數似然函數得到模型分布的參數。怎么辦呢？這就是EM算法可以派上用場的地方了。

EM算法解決這個的思路是使用：啟發式的迭代方法

• 既然我們無法直接求出模型分布參數，那么我們可以先猜想隱含數據（EM算法的E步），接著基于觀察數據和猜測的隱含數據一起來極大化對數似然，求解我們的模型參數（EM算法的M步)。
• 由于我們之前的隱含數據是猜測的，所以此時得到的模型參數一般還不是我們想要的結果。不過沒關系，我們基于當前得到的模型參數，繼續猜測隱含數據（EM算法的E步），然后繼續極大化對數似然，求解我們的模型參數（EM算法的M步)。以此類推，不斷的迭代下去，直到模型分布參數基本無變化，算法收斂，找到合適的模型參數。

從上面的描述可以看出，EM算法是迭代求解最大值的算法，同時算法在每一次迭代時分為兩步，E步和M步。一輪輪迭代更新隱含數據和模型分布參數，直到收斂，即得到我們需要的模型參數。

一個最直觀了解EM算法思路的是K-Means算法，見之前寫的K-Means聚類算法原理。在K-Means聚類時，每個聚類簇的質心是隱含數據。我們會假設K個初始化質心，即EM算法的E步；然后計算得到每個樣本最近的質心，并把樣本聚類到最近的這個質心，即EM算法的M步。重復這個E步和M步，直到質心不再變化為止，這樣就完成了K-Means聚類。

當然，K-Means算法是比較簡單的，實際中的問題往往沒有這么簡單。上面對EM算法的描述還很粗糙，我們需要用數學的語言精準描述。

EM算法的推導

Jensen不等式定義：如果f是凸函數,X是隨機變量，那么E[f(X)]≥f(E(X))。相反，如果f式凹函數，X是隨機變量，那么E[f(X)]≤f(E[X])。當且僅當X是常量時，上式取等號，其中E[X]表示X的期望。（設f是定義域為實數的函數，如果對于所有的實數X，f(X)的二次導數大于等于0，那么f是凸函數。相反，f(X)的二次導數小于0，那么f是凹函數。）

?

建議最好能夠手寫推導一遍，加深理解！

EM算法流程

EM算法的收斂性思考

EM算法的流程并不復雜，但是還有兩個問題需要我們思考：

　1） EM算法能保證收斂嗎？

　2） EM算法如果收斂，那么能保證收斂到全局最大值嗎？　　

　　首先我們來看第一個問題, EM算法的收斂性。要證明EM算法收斂，則我們需要證明我們的對數似然函數的值在迭代的過程中一直在增大。即

第二個問題：

?

從上面的推導可以看出，EM算法可以保證收斂到一個穩定點，但是卻不能保證收斂到全局的極大值點，因此它是局部最優的算法，當然，如果我們的優化目標L(θ,θj)是凸的，則EM算法可以保證收斂到全局最大值，這點和梯度下降法這樣的迭代算法相同。至此我們也回答了上面提到的第二個問題。

EM算法進一步思考

如果我們從算法思想的角度來思考EM算法，我們可以發現我們的算法里已知的是觀察數據，未知的是隱含數據和模型參數，

在E步，我們所做的事情是固定模型參數的值，優化隱含數據的分布，而在M步，我們所做的事情是固定隱含數據分布，優化模型參數的值。比較下其他的機器學習算法，其實很多算法都有類似的思想。比如SMO算法（支持向量機原理(四)SMO算法原理），坐標軸下降法(Lasso回歸算法： 坐標軸下降法與最小角回歸法小結), 都使用了類似的思想來求解問題。

EM算法優缺點

優點

• 聚類。

• 算法計算結果穩定、準確。

• EM算法自收斂，既不需要事先設定類別，也不需要數據間的兩兩比較合并等操作。

缺點

• 對初始化數據敏感。

• EM算法計算復雜，收斂較慢，不適于大規模數據集和高維數據。

• 當所要優化的函數不是凸函數時，EM算法容易給出局部最優解，而不是全局最優解。

EM算法實現

高斯混合模型(GMM)使用高斯分布作為參數模型，利用期望最大(EM)算法進行訓練。下列代碼利用高斯混合模型確定iris聚類。

import matplotlib as mpl import matplotlib.pyplot as plt import numpy as npfrom sklearn import datasets from sklearn.mixture import GaussianMixture from sklearn.model_selection import StratifiedKFoldcolors = ['navy', 'turquoise','darkorange']def make_ellipses(gmm, ax):for n, color in enumerate(colors):if gmm.covariance_type == 'full':covariances = gmm.covariances_[n][:2, :2]elif gmm.covariance_type == 'tied':covariances = gmm.covariances_[:2, :2]elif gmm.covariance_type == 'diag':covariances = np.diag(gmm.covariances_[n][:2])elif gmm.covariance_type == 'spherical':covariances = np.eye(gmm.means_.shape[1]) * gmm.covariances_[n]v, w = np.linalg.eigh(covariances)u = w[0] / np.linalg.norm(w[0])angle = np.arctan2(u[1], u[0])angle = 180 * angle / np.pi # convert to degreesv = 2.* np.sqrt(2.) * np.sqrt(v)ell = mpl.patches.Ellipse(gmm.means_[n, :2], v[0], v[1],180 + angle,color=color)ell.set_clip_box(ax.bbox)ell.set_alpha(0.5)ax.add_artist(ell)iris = datasets.load_iris()# Break up the dataset into non-overlapping training (75%) # and testing (25%) sets.skf = StratifiedKFold(n_splits=4)# Only take the first fold.train_index, test_index = next(iter(skf.split(iris.data, iris.target)))X_train = iris.data[train_index] y_train = iris.target[train_index] X_test = iris.data[test_index] y_test = iris.target[test_index]n_classes = len(np.unique(y_train))# Try GMMs using different types of covariances.estimators = dict((cov_type, GaussianMixture(n_components=n_classes,covariance_type=cov_type, max_iter=20 , random_state=0)) for cov_type in ['spherical','diag','tied','full'])n_estimators = len(estimators) plt.figure(figsize=(3* n_estimators // 2, 6)) plt.subplots_adjust(bottom=.01, top=0.95, hspace=.15, wspace=.05,left=.01, right=.99)for index, (name, estimator) in enumerate(estimators.items()):# Since we have class labels for the training data, we can # initialize the GMM parameters in a supervised manner.estimator.means_init = np.array([X_train[y_train == i].mean(axis=0)for i in range(n_classes)])# Train the other parameters using the EM algorithm.estimator.fit(X_train)h = plt.subplot(2, n_estimators // 2, index + 1)make_ellipses(estimator, h) for n, color in enumerate(colors):data = iris.data[iris.target == n]plt.scatter(data[:, 0], data[:, 1], s=0.8,color=color,label=iris.target_names[n])# Plot the test data with crossesfor n, color in enumerate(colors):data = X_test[y_test == n]plt.scatter(data[:, 0], data[:, 1], marker='x', color=color)y_train_pred = estimator.predict(X_train)train_accuracy = np.mean(y_train_pred.ravel() == y_train.ravel()) * 10plt.text(0.05, 0.9, 'Train accuracy: %.1f' % train_accuracy,transform=h.transAxes)y_test_pred = estimator.predict(X_test)test_accuracy = np.mean(y_test_pred.ravel() == y_test.ravel()) * 100plt.text(0.05, 0.8, 'Test accuracy: %.1f' % test_accuracy, transform=h.transAxes)plt.xticks(())plt.yticks(())plt.title(name)plt.legend(scatterpoints=1, loc='lower right', prop=dict(size=12)) plt.show()

在圖上，訓練數據顯示為點，而測試數據顯示為十字。?iris dataset是四維的。這里僅顯示前兩個維度，因此一些點在其他維度上分開。

參考：Sklearn官網GMM模塊http://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_covariances.html#sphx-glr-auto-examples-mixture-plot-gmm-covariances-py

? ? ? ? ? ?Pinard_EM算法原理總結

?

總結

以上是生活随笔為你收集整理的EM算法推导的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。