EM演算法簡介

首先上一段EM演算法的wiki定義:

expectation–maximization (EM) algorithm is an iterative method to find maximum likelihood（MLE） or maximum a posteriori (MAP) estimates of parameters in statistical models, where the model depends on unobserved latent variables.

就是EM演算法是: 一種迭代式的演算法，用於含有隱變數的概率引數模型的最大似然估計或極大後驗概率估計.

網上已經有很多很優秀的部落格講EM演算法的了，再次就不贅述了，只複述一些關鍵性的步驟，相關連結見本文參考部分.

(1) 起因: 給定一系列樣本，求解含有隱變數的極大似然估計(MLE)

其中z表示隱變數.
由於隱變數的存在，無法直接使用MLE去求解theta，EM的策略是先建立極大似然函式的下界(E-Step)，然後去優化下界逼近原始的極大解(M-Step)，不停迭代直到收斂到區域性最優解.

(2) 求解: 演算法推導

Qi表示隱變數z的分佈，需要滿足條件: ，比如要將班上學生聚類，假設隱藏變數z是身高，那麼Qi就是連續的高斯分佈，如果按照隱藏變數是男女，那麼就是伯努利分佈.

主要是公式2到公式3比較難懂，使用的是Jensen不等式，具體可以看這篇部落格有詳細的數學解釋，此處不贅述.

(3) 結論: 演算法總結

公式3表示是對極大似然函式求下界，此時我們假定theta已近給定，通過調整Qi的值使得下界不斷的上升去逼近真實值. 當不等式變成等式的時候表示已經調整到和真實值一樣的水平了，由Jensen不等式知道此時的隨機變數是常數C，即:

進一步推導得到:

得到第一個重要的結論:
theta已知的情況下，使得下界提升的Qi就是後驗概率，解決了Qi如何選擇的問題，其實這就是E-Step.

在找到使得下界提升的Qi之後，固定住Qi，M-Step就是使用MLE極大化此時的下界.

總結下就是:

套路就是: 首先猜下隱類別變數z，之後更新其它引數(theta)

圖解就是:

當收斂到theta*時或者||theta(t+1)-theta(t)|| < thresh就可以迭代停止了. 從演算法的過程來看，EM演算法對初始值敏感同時不能保證收斂到全域性最優解.
至於後續的證明EM演算法的收斂性，大家看我參考處的相關部落格連結或者李航博士的<<統計學習方法>>一書第9章有詳細的證明.

EM演算法的應用

GMM

GMM(Gaussian Mixture Model)就是指對樣本的概率密度（density estimation）分佈進行估計，而估計採用的模型是多個高斯模型的加權和，其中的每個高斯模型就代表了一個類(Cluster). 實際分佈中可以把模型定義為任何分佈的mixture model，為何是高斯混合模型呢? 原因如下兩點:

• 計算比較方便
• 理論任意多的高斯分佈可以近似任意概率分佈

問題簡化為:

隨機變數X是由K個高斯分佈混合而成，各個高斯分佈的權重(概率)是Φi, 各個高斯分佈的均值和方差為μi, ∑i. 觀測到隨機變數X的一系列樣本，估計引數Φ, μ, ∑.

和EM演算法之前的引入一樣，隱含類別標籤用Zi表示，表示樣本屬於類別Zi，可以假定Zi服從多項式分佈，即:

換句話來說就是第j個cluster的權重是Φj.

假設有K個類別(cluster). 假定在給定Zi後，Xi服從高斯分佈，即:

聯合概率分佈是:

故此時的極大似然函式是:

參考EM演算法的套路，首先猜測隱類別變數z，然後更新其它引數(Φ, μ, ∑).

Wji表示第i個數據點屬於第j個cluster的概率.
具體的Wji的計算可以使用貝葉斯公式:

sklearn中的GMM

The GaussianMixture object implements the expectation-maximization (EM) algorithm for fitting mixture-of-Gaussian models.

可以看出是用EM演算法求解的GMM. 官方有個示例, 示例地址是使用EM演算法來進行density estimation的.

程式碼直接貼上來，如下:

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

print(__doc__)

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()) * 100
    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()

程式碼大意是使用不同的covariance型別({‘full’, ‘tied’, ‘diag’, ‘spherical’})，來觀察GMM對iris資料集的聚類效果. iris資料集由150個樣本組成，每個樣本的特徵是4維，3個類別(setosa,versicolor,virginica).

結果如下:

EM還有用在DGM(Bayesian network)中的，這些就比較高深了，暫時還沒做了解，以後再補.