EM算法推导

- 日理万妓 2022-05-18 07:37 233阅读 0赞

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算法的推导 ###

![70][]

**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是凹函数。）

![70 1][]

建议最好能够手写推导一遍，加深理解！

### EM算法流程 ###

![70 2][]

### EM算法的收敛性思考 ###

EM算法的流程并不复杂，但是还有两个问题需要我们思考：

1） EM算法能保证收敛吗？

2） EM算法如果收敛，那么能保证收敛到全局最大值吗？

首先我们来看第一个问题, EM算法的收敛性。要证明EM算法收敛，则我们需要证明我们的**对数似然函数的值在迭代的过程中一直在增大**。即

![70 3][]

第二个问题：

![70 4][]

从上面的推导可以看出，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()

![70 5][]

在图上，训练数据显示为点，而测试数据显示为十字。 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][http_scikit-learn.org_stable_auto_examples_mixture_plot_gmm_covariances.html_sphx-glr-auto-examples-mixture-plot-gmm-covariances-py]

Pinard\_EM算法原理总结

[70]: /images/20220518/d675cdb1d6bc41b5bf30ac9a5837275f.png
[70 1]: /images/20220518/7a1978f2c3f647b294dde0845cd04b07.png
[70 2]: /images/20220518/436ad67d7c274b3b926b4408181f8682.png
[70 3]: /images/20220518/e693c82da8bb47afae4e17119426f2cf.png
[70 4]: /images/20220518/9f9731679b9549d5b08231bb367e93cb.png
[70 5]: /images/20220518/fa1e2f4f0c3f4849b2ae76b35fb959a5.png
[http_scikit-learn.org_stable_auto_examples_mixture_plot_gmm_covariances.html_sphx-glr-auto-examples-mixture-plot-gmm-covariances-py]: http://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_covariances.html#sphx-glr-auto-examples-mixture-plot-gmm-covariances-py