什么是高斯混合模型

作者：Oscar Contreras Carrasco

翻译：老齐

图书推荐：《数据准备和特征工程》

机器学习可以分为两个主要领域：有监督学习和无监督学习。两者的主要区别在于数据的性质以及处理数据的方法。聚类是一个无监督学习的算法，利用这个算法可以从数据集里找到具有共性的点簇。假设我们有一个如下所示的数据集：

我们的工作是找到看起来很接近的点簇(聚类)。按照图中所示的分布，可以清楚地识别出两个聚类，我们将其分别用蓝色和红色标出:

请注意，我们现在引入了一些附加的符号。这里，μ1和μ2是每个聚类的质心，也是识别每个聚类的参数。一种流行的聚类算法被称为K-means(K均值)，它用遵循迭代的方法来更新每个聚类的参数。更具体地说，它要做的是计算每个聚类的平均值（或质心），然后计算质心到每个数据点的距离，后者被标记为聚类的一部分，这个聚类是由其最近的质心来标识的。这个过程会重复，直到满足某些收敛条件。例如，当我们看到聚类的赋值没有进一步的变化时。

K-means(K均值)的一个重要特点是它是一种硬聚类方法，它将每个点与一个(且仅与一个)聚类相关联。这种方法的一个局限性是没有不确定性度量标准或概率来告诉我们一个数据点与一个特定的聚类的关联程度。那么，如果使用软聚类而不是硬聚类，效果会怎么样呢？这正是高斯混合模型（简称GMMs）所要尝试的。现在我们来进一步探讨这个方法。

定义

高斯混合模型由多个高斯函数组成，每个高斯甘薯由

标识，其中

是数据集的簇数(聚类数)。混合模型中的每个高斯

由以下参数组成：

• 定义其中心的平均值 μ。
• 定义其宽度的协方差 ∑。这相当于多变量场景中椭球体的维数。
• 混合概率 π，定义高斯函数的大小。

现在，我们用图形说明这些参数：

在这里，我们可以看到三个高斯函数，因此

。每个高斯函数分别解释了三个可用聚类中包含的数据。混合系数本身就是概率，必须满足以下条件：

如何确定这些参数的最佳值呢？为了实现这一目标，必须确保每个高斯函数所对应的数据点都属于对应的一个聚类，这正是最大似然法的作用。

一般来说，高斯密度函数由以下公式给出：

其中x代表数据点，D是每个数据点的维数。μ和∑分别是均值和协方差。如果我们有一个由N=1000个三维点组成的数据集（D=3），那么x就是1000×3矩阵。μ是1×3矢量，∑是3×3矩阵。为了便于之后的使用，我们还会发现该方程的对数是有用的，其公式如下：

如果我们根据平均值和协方差对该方程进行微分，然后使其等于零，就能够找到这些参数的最佳值，并且解决方案将对应于该设置的最大似然估计（MLE）。但是，因为我们处理的高斯函数不仅仅是一个，而是多个，当我们找到整个混合的参数时，情况会变得有点复杂。因此，需要介绍其他方面的一些知识，这也是下一节中将要讨论的。

导数

我们现在要介绍一些附加的符号。警告一句：数学来了！别担心。为了更好地理解推导过程，我将尽量保持符号的清晰。首先，假设我们想知道数据点

来自高斯分布

的概率是多少，可以将其表示为：

它的意思是：“给定一个数据点x，它来自高斯分布 k 的概率是多少？” 在本例中，z是一个潜在变量，它只接受两个可能的值。当x来自高斯k时，z的值为1，否则z的值为0。实际上，我们在现实中并没有看到这个z变量。但是，了解数据点x来自高斯k的概率将有助于我们确定高斯混合参数，正如我们后面将会讨论的那样。

同样，我们可以说：

这意味着：观察到来自高斯k的点的总概率实际上等于高斯的混合系数。这是有道理的，因为高斯分布越大，我们期望的概率就越高。现在让z成为所有可能的潜在变量z的集合，因此：

我们预先知道：每个z独立于其他z出现，并且当k等于k点所在的聚类时，它们只能取1的值。因此：

现在，如果我们的数据来自高斯k，那么如何找到观测数据的概率呢？原来它实际上是高斯函数本身！按照定义p(z)的相同逻辑，可以说：

现在你可能会问，为什么要这么做？还记得最初的目标吗？它就是：已知观测值为x, 确定z的概率。结果证明，刚刚推导的方程，连同贝叶斯规则，将帮助我们确定这个概率。由概率的乘积法则可知：

看来我们现在有进展了，等号右边就是最新的结果，也许你已经预料到：我们将使用贝叶斯方法来获取最终需要的概率。但是，首先需要

，而不是

。怎么才能去掉z呢？是的，你猜对了。边缘化！我们只需要把z上的项求和，因此

这是定义高斯混合模型的方程，可以清楚地看到它依赖于在前面提到的所有参数！为了确定最佳值，需要确定模型的最大似然。可以找到所有观测值

的联合概率，定义如下：

就像我们对原始高斯密度函数所做的那样，我们将对数应用到方程的两侧：

太好了！为了找到高斯混合模型的最佳参数，我们所要做的就是把这个方程和参数相比较，任务完成了，对吧？等一等！没有这么快。还有一个问题。我们可以看到，有一个对数影响了第二次求和。计算这个表达式的导数，然后求解参数，这是非常困难的！

怎么办？需要用迭代的方法来估计参数。还记得在已知x的情况下，如何找出z的概率吗？行动起来吧，因为在这一点上，我们已经有了定义这个概率的所有条件。

根据贝叶斯定律，得：

从之前的推导中，我们知道:

现在用前面的等式来替换它们：

这就是我们一直在寻找的！继续下去，会经常看到这种表达方式。接下来，将继续讨论一种方法，它将帮助我们很容易地确定高斯混合模型的参数。

期望最大化算法

现在，已经推导出了一些概率表达式，我们将发现这些表达式在确定模型参数时很有用。然而，在上一节中，仅仅通过计算（3）来找到这样的参数是非常困难的。幸运的是，有一种迭代方法可以用来达到这个目的。它被称为期望最大化，或简称EM算法。它被广泛应用于目标函数具有复杂性的优化问题，例如我们刚刚在GMM案例中遇到的问题。

设模型参数为

现在我们来定义一般EM算法将要遵循的步骤。

步骤1:相应地初始化 θ 。例如，我们可以使用之前运行k均值获得的结果，作为算法的一个良好起点。

步骤2(期望步骤):评估

实际上我们已经求出了p(z|x, θ)。还记得我们在上一节中使用的γ 表达式吗? 为了表达得更清晰，我们把之前的方程(4)带入这里:

对于高斯混合模型，期望步骤归结为利用旧的参数值计算(4)中γ 的值。现在将(4)替换为(5)，得到:

听起来不错，但是我们仍然缺少p(X, Z|θ* )。我们怎样才能找到它呢? 实际上这并不难。它正好是模型的完全似然性，包括X和Z。我们可以用如下表达式来求它:

它是计算所有观测值和潜在变量的联合概率的结果，是对p(x)的初始推导的扩展。这个表达式的对数是：

很好！终于摆脱了这个麻烦的对数，它影响了(3)中的求和。现在万事俱备了。对我们来说，更容易估算参数的办法是：使Q相对于参数最大化。但是我们将在最大化步骤中处理这个问题。此外，请记住，每次求和时，潜在变量z仅为1。有了这些知识，我们就可以很容易地在推导过程中消除它。

最后，我们把(7)替换成(6)，得到：

在最大化步骤中，我们将得到修正后的混合参数。为此，需要使Q成为一个受限最大化问题，因此将向(8)添加一个拉格朗日乘子。现在我们回顾一下最大化步骤。

步骤3(最大化步骤): 使用下列等式得到修正的参数 θ*

其中

这就是我们在上一步中得到的结果。然而，Q也应该考虑到限制因素：所有的π值之和应为1。为此，需要添加一个合适的拉格朗日乘数。因此，我们应该将(8)改写为:

现在我们可以用极大似然很容易地确定参数。取Q对π的导数，设它为零

然后，通过重新排列这些项，并将k上的求和应用于方程的两侧，我们得到：

由(1)可知，所有混合系数π 的总和等于1。另外，我们知道，把k上的概率γ 加起来也会得到1。利用这个结果，我们可以求解出π:

同样，如果我们将Q与μ和∑相区别，让导数等于零，然后利用我们定义的对数似然方程（2）求解参数，我们得到：

就是这样！然后，我们将使用这些修正值来确定下一个EM迭代中的γ，以此类推，直到我们看到似然值的一些收敛性。我们可以用方程（3）来监测每一步的对数似然，并且我们总是保证达到局部极大值。

如果我们能用编程语言来实现这个算法，那就太好了，不是吗？接下来，我们将看到我所提供的Jupyter笔记本的部分内容，这样你就可以看到GMMs在Python中的工作实现。

在Python中的实现

补充说明：在https://bit.ly/2MpiZp4可以看到作为Jupyter笔记本的完整的实现。

我在这个练习中使用了Iris数据集，主要是为了简单和快速的训练。在之前的推导中已知：EM算法遵循迭代的方法来寻找高斯混合模型的参数。我们的第一步是初始化参数。在这种情况下，可以使用K均值的值来满足这个目的。相关的Python代码如下所示：

def initialize_clusters(X, n_clusters):
    clusters = []
    idx = np.arange(X.shape[0])
    
    kmeans = KMeans().fit(X)
    mu_k = kmeans.cluster_centers_
    
    for i in range(n_clusters):
        clusters.append({
            'pi_k': 1.0 / n_clusters,
            'mu_k': mu_k[i],
            'cov_k': np.identity(X.shape[1], dtype=np.float64)
        })
        
    return clusters

接下来，执行期望步骤：

相应的Python代码如下所示：

def expectation_step(X, clusters):
    totals = np.zeros((X.shape[0], 1), dtype=np.float64)
    
    for cluster in clusters:
        pi_k = cluster['pi_k']
        mu_k = cluster['mu_k']
        cov_k = cluster['cov_k']
        
        gamma_nk = (pi_k * gaussian(X, mu_k, cov_k)).astype(np.float64)
        
        for i in range(X.shape[0]):
            totals[i] += gamma_nk[i]
        
        cluster['gamma_nk'] = gamma_nk
        cluster['totals'] = totals
    
    for cluster in clusters:
        cluster['gamma_nk'] /= cluster['totals']

注意，为了计算求和，只需要利用分子中的项并相应地做除法。

然后，我们得到最大化步骤，在这里我们计算：

相应的Python代码如下：

def maximization_step(X, clusters):
    N = float(X.shape[0])
  
    for cluster in clusters:
        gamma_nk = cluster['gamma_nk']
        cov_k = np.zeros((X.shape[1], X.shape[1]))
        
        N_k = np.sum(gamma_nk, axis=0)
        
        pi_k = N_k / N
        mu_k = np.sum(gamma_nk * X, axis=0) / N_k
        
        for j in range(X.shape[0]):
            diff = (X[j] - mu_k).reshape(-1, 1)
            cov_k += gamma_nk[j] * np.dot(diff, diff.T)
            
        cov_k /= N_k
        
        cluster['pi_k'] = pi_k
        cluster['mu_k'] = mu_k
        cluster['cov_k'] = cov_k

请注意，为了简化计算，我们使用了：

最后，我们还进行了对数似然计算：

相应的Python代码是这样的：

def get_likelihood(X, clusters):
    likelihood = []
    sample_likelihoods = np.log(np.array([cluster['totals'] for cluster in clusters]))
    return np.sum(sample_likelihoods), sample_likelihoods

已经在期望步骤中预先计算了第二个求和的值，所以我们在这里就利用它。此外，创建图表来查看似然性的进展情况总是很有用的。

可以清楚地看到算法在大约20次之后收敛。EM保证在给定的过程迭代次数后将达到局部最大值。

最后，作为实现的一部分，我们还生成一个动画，向我们展示每次迭代后聚类设置是如何改进的。

注意GMM如何改进质心，这些质心是通过K均值估计的。当我们收敛时，每个聚类的参数值不会进一步改变。

高斯混合模型是一种非常强大的工具，广泛应用于涉及数据聚类的各种任务中。

原文链接：https://towardsdatascience.com/gaussian-mixture-models-explained-6986aaf5a95