EM算法推導

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 np


from  sklearn import datasets
from  sklearn.mixture import GaussianMixture
from  sklearn.model_selection import StratifiedKFold


colors = ['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 degrees
      v = 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 crosses
 
  for 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()) * 10
        plt.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()) * 100

        plt.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算法原理總結