机器学习

GMM模型原理与实践

1. 模型

我们假设这篇文章的读者对于EM算法有着一定了解。

1.1 原理

在之前的文章EM算法原理与推导中我们提到：

假设我们有数据集\(X={x_{(1)},x_{(2)},…,x_{(m)}}\)，我们想利用生成式模型建模\(P(x;\theta)\)。…… 认为数据来自于多个分布，只不过我们无法观测到数据被哪一个子分布生成，因此不能直接利用公式直接计算解析解。EM算法就是用来处理这种具有“隐标签”的问题的。

同样，在PRML一书第九章第一句话也给出了提纲挈领的表述：

如果我们定义观测变量和潜在变量的一个联合概率分布，那么对应的观测变量本身的概率分布可以通过求边缘概率的方法得到。这使得观测变量上的复杂的边缘概率分布可以通过观测变量与潜在变量组成的扩展空间上的更加便于计算的联合概率分布来表示。因此，潜在变量的引入使得复杂的概率分布可以由简单的分量组成。

高斯混合模型就是这种情况的特例。我们以比较简单的语言整理一下这两段话。（不要在意把我自己的文章中的话和机器学习圣经中的话放到一起这件事…）
关键是弄明白这三个粗体的“概率”的意思，我们是在为无标签数据建立无监督的生成式模型\(P(x;\theta)\)。在这里，所谓模型，就是概率。$$$ P(x;\theta)=\sum_zP(x^{(i)},z^{(i)};\theta) $$$，就是指“观测变量本身的概率分布可以通过求边缘概率的方法得到。”，用模型的方式说就是混合模型的就是多个模型的累加和。而边缘概率又是对联合概率中某一随机变量的求和，之所以“定义观测变量和潜在变量的一个联合概率分布”，是因为联合概率分布就是子模型，而子模型好求。而单个父模型扩展成多个子模型的原因是，我们增加了一个变量，这个变量是隐变量，从而人为增加了一个变量维度。

2. 参数估计

子模型（联合分布）是高斯分布：

$$ P(X, Y=i)=P(X|Y=i) \times P(Y=i) $$

$$ =\frac{1}{(2\pi)^{n/2}|\Sigma_i|^{1/2}}exp\Big(-\frac{1}{2}(x-\mu_i)^T\Sigma_i^{-1}(x-\mu_i)\Big) \times \phi_i $$

总模型（边缘分布）是子模型在隐变量上的累加：
$$ P(X)=\sum_iP(X, Y=i) $$

既然我们增加了隐变量，极大似然就不能直接用了，然后就能套到EM算法中了。

批改于(2018-06-01)：上面这句话太草率了，其实不是MLE不能用了，直接对混合模型的对数似然求梯度使其等于零依然可以计算，只不过计算得到的参数估计值不再是一个常规的解析式了，而是依赖于对隐变量的后验概率，详见EM算法原理与推导(续)。

$$ l(\theta)=\sum_{i=1}^mlogP(x^{(i)};\theta)\qquad$$
$$ l(\theta)=\sum_{i=1}^mlog\sum_zP(x^{(i)},z^{(i)};\theta)\quad$$

根据之前EM算法中进行到最后的推导，E步需要求隐标签的后验分布，这个是好求的。我们已知\(P(X|y)\)的正态分布计算方式，那么通过贝叶斯公式就可以求得反向的条件概率。用$$$\omega$$$表示我们计算得到的隐变量的后验概率，即$$$Q$$$，则

$$ \omega_j^{(i)}=Q_i(z^{(i)}=j)=P(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma) $$

然而在EM算法原理与推导中，M步的时候，留下了一个疑问：M步的计算公式好算吗，为什么好算？

$$ \theta := argmax_\theta J(Q, \theta) $$

2.2 参数估计：$$$\mu$$$(M-Step)

至少目前看来，在高斯模型中是这样的：直接累加的模型中有隐变量我们不能求解MLE，如果使用了EM算法，真正M步可以计算的原因不是公式的形式，而是我们在E步中确定下来的隐变量的后验概率，那么我们就可以把这个值带入进去进行计算了。联合分布的概率公式我们也明确的知道。即：

$$ J(Q, \theta)=\sum_i^m\sum_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} $$
$$ J(Q, \theta)=\sum_{i=1}^m\sum_{j=1}^kQ_i(z^{(i)}=j)log\frac{P(x^{(i)}|z^{(i)}=j;\mu,\Sigma)P(z^{(i)}=j;\phi)}{Q_i(z^{(i)})} $$
$$ J(Q, \theta)=\sum_{i=1}^m\sum_{j=1}^k\omega_j^{(i)}log\frac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp\Big(-\frac{1}{2}(x-\mu_j)^T\Sigma_j^{-1}(x-\mu_j)\Big)*\phi_j}{\omega_j^{(i)}} \qquad(1)$$

令$$$J(Q, \theta)$$$对$$$\mu_l$$$的梯度为零：
$$ \nabla_{\mu_l}\sum_{i=1}^m\sum_{j=1}^k\omega_j^{(i)}log\frac{\frac{1}{(2\pi)^{n/2}|\Sigma_j|^{1/2}}exp\Big(-\frac{1}{2}(x-\mu_j)^T\Sigma_j^{-1}(x-\mu_j)\Big)*\phi_j}{\omega_j^{(i)}} $$
$$ =-\frac12\nabla_{\mu_l}\sum_{i=1}^m\sum_{j=1}^k\omega_j^{(i)}(x-\mu_j)^T\Sigma_j^{-1}(x-\mu_j) $$
$$ =-\frac12\sum_{i=1}^m\nabla_{\mu_l}\omega_l^{(i)}(x-\mu_l)^T\Sigma_l^{-1}(x-\mu_l) $$
$$ =\sum_{i=1}^m\omega_l^{(i)}\Sigma_l^{-1}(x-\mu_l)=\sum_{i=1}^m\omega_l^{(i)}(x-\mu_l)=0$$
所以：
$$ \mu_l:=\frac{\sum_{i=1}^m\omega_l^{(i)}x^{(i)}}{\sum_{i=1}^m\omega_l^{(i)}} $$

2.3 参数估计：$$$\phi$$$(M-Step)

将（1）式中的log展开，只保留与$$$\phi$$$有关的项，则我们实际需要优化的项是：
$$ \sum_{i=1}^m\sum_{j=1}^k\omega_j^{(i)}log\phi_j $$

优化这一项有一点需要注意，$$$log\phi$$$是一个单调函数，不能直接MLE，而限制它的是另外的条件：$$$\phi$$$是一个分布，要求和为一。因此带约束的优化问题，考虑使用拉格朗日乘子法：

$$ \mathcal{L}(\phi)=\sum_{i=1}^m\sum_{j=1}^k\omega_j^{(i)}log\phi_j+\beta(\sum_{j=1}^k\phi_j-1)$$
$$ \frac{\partial}{\partial\phi_l}=\frac1{\phi_l}\sum_{i=1}^m\omega_l^{(i)}+\beta=0 $$
$$ \therefore \qquad \phi_l:=\frac{\sum_{i=1}^m\omega_l^{(i)}}{-\beta} \qquad$$

$$$\beta$$$的求法这里不多说，是Lagrange Multiplier的基本方法，就是把上面的$$$\phi_l$$$的表达式带入约束条件即求和为一，可以解得$$$-\beta=m$$$，因此：
$$ \phi_l:=\frac{\sum_{i=1}^m\omega_l^{(i)}}{m} $$

注意：其实还有$$$\phi$$$必须大于零的要求，我们没有考虑是因为这里即使只考虑等式约束，解得的解也会自动满足不等式约束，深入的原因在于广义拉格朗日乘子法。

2.4 参数估计：$$$\Sigma$$$(M-Step)

协方差矩阵的推导不多做说明，可以参考GDA原理与推导，过程类似，这里只给出解析式。

$$ \mu_l:=\frac{\sum_{i=1}^m\omega_l^{(i)}(x^{(i)}-\mu_l)(x^{(i)}-\mu_l)^T}{\sum_{i=1}^m\omega_l^{(i)}} $$

sklearn实践

sklearn.mixture中实现了GMM和变分贝叶斯GMM。主要实践一下GMM，下面是从sklearn官方示例演变的example。详见参考内容1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import matplotlib  as mpl
from matplotlib import  pyplot as plt
import numpy as np
from sklearn.datasets import make_classification
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)

X, y = make_classification(n_samples=600, n_features=2, n_redundant=0, n_informative=2, n_classes=3, n_clusters_per_class=1, random_state=15, flip_y=0)


# Break up the dataset into non-overlapping training (75%) and testing
# (25%) sets.
skf = StratifiedKFold(n_splits=5)
# Only take the first fold.
train_index, test_index = next(iter(skf.split(X, y)))
X_train = X[train_index]
y_train = y[train_index]
X_test = X[test_index]
y_test = y[test_index]

def show_raw_data():
	# 以散点图的形式展示原始数据
    plt.scatter(X_train[:, 0], X_train[:, 1], marker='.', c=y_train, s=10)
    plt.scatter(X_test[:, 0], X_test[:, 1], marker='o', c=y_test, s=25)
    plt.show()

# show_raw_data()

# 注意这个n_classes也是把GMM当做监督学习来用，通过标签确定组分个数
# 然而作为非监督学习的GMM，这个参数的确定才是困难
n_classes=len(np.unique(y))
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.

    # 这里是关键！将无监督算法当做有监督算法的话，用带标签数据去初始化参数
    # 如果把这个去掉，不要改变生成的数据（包括n_sample）和gmm模型的初始化（主要是random_state， 会发现结果完全改变，原因想一下就能明白）
    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 = X_train[y_train == n]
        plt.scatter(data[:, 0], data[:, 1], s=0.8, color=color)

    # 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()

def show_predict_area():
    x_grid = np.linspace(-5,5,101)
    y_grid = np.linspace(-5,5,101)
    x_grid, y_grid = np.meshgrid(x_grid, y_grid)
    pred_area = estimators['full'].predict(tuple(zip(x_grid.flatten(),y_grid.flatten())))

    plt.scatter(x_grid.flatten(), y_grid.flatten(), c=tuple(map(lambda x:colors[x], pred_area)), alpha=0.04)
    plt.scatter(X_train[:, 0], X_train[:, 1], marker='.', c=tuple(map(lambda x:colors[x], y_train)), s=10)
    plt.scatter(X_test[:, 0], X_test[:, 1], marker='o', c=tuple(map(lambda x:colors[x], y_test)), s=25)
    plt.show()

# show_predict_area()

3. 补充

3.1 过拟合问题

GMM可能会导致奇异解，即一个component只拟合了一个数据点，在最大似然的过程中，就会导致该组分的方差趋近于无穷小，这其实就是发生了过拟合问题。如果我们使用贝叶斯方法，这种情况就不会出现，即为各组分添加先验概率，然后利用EM算法进行极大后验估计。
另外的启发式方法是如果检测到高斯组分收缩到一个点，就重启该组分，为其参数随机赋值。

3.2 可区分问题（identifiability）

由于是无监督问题，K个分量的混合模型有K!个等价的解，HMM模型中也有该问题，注意下面的代码中是如何用监督方式的初始化来避免该问题的。

4. 对比K-Means算法

与Kmeans算法相比，GMM在收敛前迭代次数更多，计算量也更大，通常可以使用Kmeans找到一个合适的聚类结果来初始化GMM。
KMeans算法和GMM模型有着很强的关系，如果认为所有的高斯成分共享协方差矩阵$$$\epsilon I$$$，那么reponsibility就可以简化为:
$$\frac{\phi_kexp\Big(-\frac{1}{2\epsilon}||x-\mu_k||^2\Big) }{\sum_j\phi_jexp\Big(-\frac{1}{2\epsilon}||x-\mu_j||^2\Big) }$$
上式中，如果$$$\epsilon$$$趋近于0，可以想见，中心离该数据点最近的高斯组分的后验概率会趋近于1，其他的则会趋近于0，而得到对数据的“硬分配”，而最大化对数似然也会转变为最小化均方误差（因为共享协方差）。

参考：

  1. sklearn 文档：GMM covariances
本站总访问量次