高斯混合模型（GMM）与期望最大化算法（EM）

峡谷相对论

已于 2023-05-06 11:16:25 修改

阅读量5.3k

点赞数 11

分类专栏：算法文章标签：概率论机器学习最大期望算法语音识别

于 2020-06-25 23:00:13 首次发布

本文链接：https://blog.csdn.net/qq_33552519/article/details/106963417

版权

算法专栏收录该内容

19 篇文章 22 订阅

订阅专栏

1.1 极大似然估计

一般来说，多个相互独立的事件都发生了的概率是非常低的。当概率分布类型已知而参数待定时，我们更加倾向于认为最合适的参数应该会使得这个组合事件所发生的概率在所有可选的参数中达到最大。假设我们有一个随机信源 ${S}$ ，其分布类型已知，但参数 ${\bf{\theta }}$ 未知。那么，如果我们对信源 ${S}$ 进行多次独立同分布采样，并且得到了一系列的随机样本，其样本值为 ${\bf{X}} = \left\{ {{{\bf{x}}_1},{\rm{ }}{{\bf{x}}_2},{\rm{ }} \ldots ,{\rm{ }}{{\bf{x}}_N}} \right\}$ ，我们就有理由认为最优的分布参数会使得样本集 ${\bf{X}}$ 出现的概率最大。这时，我们只需要把分布参数作为自变量，以该参数下样本集出现的概率作为因变量，进行最大化优化即可。[注：一般概率论书籍里以大写字母为随机变量，小写字母为随机变量的值，这里为了简便，以大写字母作为采样值的集合。采样值一般为向量，有时为了方便推导可能只用到标量]。因此，参数 ${\bf{\theta }}$ 的优化就可表示为式（1-1），即
1-1
这就是极大似然估计（Maximum Likelihood Estimation, MLE）算法的基本思路。因为连乘不方便求导，而且概率的连乘容易导致小数位溢出，我们令 1-2
因为 ${\log}$ 函数是严格单调的，所以两者的极值点也是一样的。那么求解式（1-3）即可得到最优的参数，
1-3
当然这通常只能保证是局部最优的。

以一元单高斯分布为例，
1-4
根据式（1-3）可得
1-5
可以看到，对于单高斯这种分布函数，通过 MLE 来求解最优参数的计算过程还是比较简单的。

1.2 高斯混合模型

上面我们讨论了 MLE 求解单高斯分布参数的方法，然而，考虑这么一种情况，我们先从一个有限整数随机数发生器中随机抽取一个数 ${y}$ ，例如 ${1, 2, ..., M}$ ，然后根据抽取到的结果决定使用具有不同参数的分布，例如高斯分布，来生成最终的输出 ${\bf{x}}$ ，其中
1-6
那么输出 ${\bf{x}}$ 的概率或概率密度可表示为
1-7
这就是所谓的多模态（Multimodal）概率模型。如果所用到的分布为高斯分布，则称为高斯混合模型（Gaussian Mixture Model, GMM），如图 1-1 所示。可以看到，如果只用一个高斯分布来描述则会造成比较大的误差。

图1-1

图1-1 高斯混合模型

这时，如果我们想通过极大似然估计来优化各个分布函数的参数，根据式（1-2）可得，
1-8
因为 ${\log}$ 函数里面包含了加法，一般来说我们是很难得到式（1-3）的解析解的，所以得借助一些非线性优化方法，例如梯度下降（上升）法等等，但这些方法还是需要用到式（1-8）的梯度，牛顿法甚至还需要用到二阶偏导数，除非是计算机自动微分，依靠人工推导还是十分麻烦的。而且，这些方法的收敛精度和收敛速度十分依赖于泰勒公式展开的近似程度以及学习速率的设置，因此具有一定的不稳定性。

然而，高斯分布有着其独特的性质。对于单高斯模型，根据 MLE 我们可以得到其两个最优估计的参数刚好为采样数据的均值和方差，即如式（1-5）所示。那么对于 GMM 来说，如果我们知道每个样本的类别，那是否也可以根据式（1-5）对每个高斯分布的参数进行估计？当然，我们是不可能知道这些样本的类别的，但是，如果每个高斯分布的参数已知，通过贝叶斯定理我们可以得到这些样本属于某个类别的后验概率，即
1-9
那么，我们就可以对各个高斯分布的参数进行软估计了，即
1-10
慢着！这绕了一圈下来不就成了鸡生蛋、蛋生鸡的问题了？在求某个样本属于各个类别的后验概率时我们需要先知道每个高斯分布的参数，但这些参数又是基于前面求得的后验概率的。因此，很自然地，我们想到了一种迭代优化的方法，即先随机选择各个高斯分布的参数，然后通过式（1-9）和（1-10）得到新的参数，依此重复，最后或许我们能够得到足够好的估计结果。注意，是“或许”！我们怎么能够保证这种迭代方法一定能够使得参数往极大似然的方向前进呢？我们这里甚至连梯度都没有提到！为了解决这个问题，我们需要引入期望最大化（Expectation Maximization, EM）算法。通过 EM 算法，我们将证明，前面所述的迭代方法是可以保证达到一个局部最优点的。这也是为什么 GMM 应用能够这么广泛的一个原因，除了这种分布在自然中广泛存在以外，其参数优化的过程也是比较简单的。

1.3 期望最大化算法

从直觉上，我们认为可通过式（1-9）和（1-10）的迭代来得到高斯混合模型参数的一个最优估计，但是我们目前没办法从理论上去证明这种迭代的有效性，而 EM 算法就是解决这个问题的法宝。对于大多数多模态系统来说，我们只能获取最终的输出 ${\bf{x}}$ ，而无法知道该输出所属的类别 ${y}$ 。因此，我们一般将 ${y}$ 称为隐含变量， ${\bf{x}}$ 则称为不完整数据（Incomplete Data），而 ${{\bf{z}} =({\bf{x}}, y)}$ 才称为完整数据（Complete Data）。EM 算法就是要解决这种包含不完整数据的极大似然估计问题。实际上，隐含变量 ${y}$ 不一定是一个多模态分布中的类别信息，其甚至可以只是某个样本值 ${\bf{x}}$ 中的某个维度或特征缺失的值，但这里我们只讨论 ${y}$ 属于类别信息的情况。

1.3.1 迭代策略

EM 算法首先定义了一种能保证收敛的迭代过程，如图 1-2 所示。给定一个需要寻找极大值点的函数 ${f(x)}$ ，找到一个辅助函数 ${A(x, x^t)}$ ，其中 ${x^t}$ 为 ${t}$ 时刻所找到的最优参数值，也就是辅助函数是与当前最优解相关的，并且满足
1-11
那么，找到 ${A(x, x^t)}$ 的极大值点，则有
1-12
这样，我们就可以保证整个迭代过程是不会恶化的，最终迭代停滞的时候就说明 ${f(x)}$ 到达了一个极大值点。要证明这一点，令
1-13
那么 ${x^t}$ 一定是 ${h(x)}$ 的一个极小值点，所以 ${f(x)}$ 和 ${A(x, x^t)}$ 在 ${x^t}$ 处的导数是相等的，也就是两函数在 ${x^t}$ 处相切。因此，迭代停滞则意味着
1-14
因为一般的优化方法只能保证 ${x^{t+1}}$ 是一个局部极大值，因此 EM 算法的结果也只能最多保证是一个局部最优解。为了获取更好的结果，通常会选择不同的初始点进行迭代，加大了找到全局最优解的几率。

图1-2

图1-2 EM 算法的迭代过程

然而怎样找到一个合适的辅助函数 ${A(x, x^t)}$ ？这似乎毫无头绪。而且我们不能让 ${A(x, x^t)}$ 的形式过于复杂，例如在 ${\log}$ 函数中包含了加法项，因为我们还是得求解 ${A(x, x^t)}$ 的极值点。当然，即便没有找到 ${A(x, x^t)}$ 的极值点，只要满足 ${A(x^{t+1}, x^{t})>A(x^t, x^t)}$ ，我们至少还是能够往好的方向前进的。这样问题最终还是落在了如何找到合适的 ${A(x, x^t)}$ 上，我们希望其能够像单模态分布的似然函数那样 ${\log}$ 函数里面不包含加法项，这样求解极值点也会简便很多。

1.3.2 Jensen不等式

在推导辅助函数 ${A(x, x^t)}$ 的具体形式之前，我们有必要先介绍一下 Jensen 不等式。首先，一个函数 ${f(x)}$ 在区间 ${[a, b]}$ 内称为凸（Convex）函数，如果
1-15
如图 1-3 所示。注意函数的凹凸从直观上是以函数以上部分的空间形状来定义的，所以数学上的凸函数看上去实际上凹形的。

图1-3

图1-3 凸函数

那么，对于一个凸函数，根据 Jensen 不等式有
1-16
我们可以使用归纳法来证明。 ${n=1}$ 时很明显成立，对于 ${n=2}$ 以上的，有
1-17
得证。因为 ${-\log(x)}$ 在 $\infty)}$ 是严格的凸函数，根据 Jensen 不等式有
1-18
利用这个不等式，我们将找到我们所需要的那个辅助函数 ${A(x, x^t)}$ 。

1.3.3 EM算法推导

回到多模态分布的参数优化，有了 Jensen 不等式和式（1-18），我们很自然地可能想改写式（1-7），即
1-19
然而这样我们无法建立起参数的迭代过程，所以还得另辟蹊径。为了构造辅助函数，我们引入另外一个概率分布 ${q(y)}$ ，其中 ${y=1, 2, ..., M}$ 。那么有
1-20
这似乎把问题复杂化了，而且看上去我们同样也没法进行参数的迭代更新。但是，因为 ${q(y)}$ 是可以自由选择的，我们或许可以通过 ${q(y)}$ 建立一个从旧参数到新参数的桥梁，从而使得迭代过程能够运行起来。

回顾式（1-11），辅助函数除了需要满足不大于原函数的条件以外，还必须在旧参数处等于原函数。因此，我们首先看看式（1-20）两侧的差距，有
1-21
也就是说，式（1-20）中的不等式两侧的差距刚好为 ${q(y)}$ 和 ${p(y|\bf{x, \Theta})}$ 的 KL 散度。我们同样可以利用 Jensen 不等式证明 KL 散度是非负的，且可以证明当且仅当 ${q(y)}$ 和 ${p(y|\bf{x, \Theta})}$ 相同的时候为 0，我们这里取其结论而略去相关证明。因为我们只要求原函数和辅助函数在旧参数处的值相等，所以我们令
1-22
于是我们得到
1-23
并且有
1-24
因此，下一步的重点是求解辅助函数的极值点，我们希望辅助函数能够让这个计算过程变得更加简单。由式（1-23）可得
1-25
那么，我们实际只需要求解 ${Q({\bf{\Theta, \Theta}}^{old})}$ 的极值点即可。而且这时 ${\log}$ 函数里面已经没有加法项了，对于高斯分布来说这能够极大地简化极值点的求解。

另外，细心的你可能已经发现了，单讨论一个样本， ${Q({\bf{\Theta, \Theta}}^{old})}$ 实际上就是 ${{\bf{z}}_i}$ 也就是 完整数据的 ${\log}$ 似然的条件期望值，英文表述为 Conditional expectation of the complete data log lokelihood，这就是 EM 算法的核心意义，也就是说将不完整数据的极大似然估计转换为完整数据的条件期望对数似然的最优化。因为完整数据的概率一般是条件概率的连乘，而不是不完整数据的在各种可能事件下的概率叠加，从而很好地避免了 ${\log}$ 函数内包含加法的情况。 ${Q({\bf{\Theta, \Theta}}^{old})}$ 优化的具体过程我们先留到后面与 GMM 的参数估计一起详细讨论，现在我们可以总结 EM 算法的基本流程如下：

(1) 选择一个初始的参数集 ${{\bf{\Theta}}^{old}}$ ，但为了避免遇到较差的局部最优解，一般会选择多组的初始参数集独立优化。
(2) E Step：根据式（1-9）计算后验概率 ${p(y|{\bf{x}}_i, {\bf{\Theta}}^{old})}$ ，构造 ${Q({\bf{\Theta, \Theta}}^{old})}$ 。
(3) M Step：求解 ${{\bf{\Theta}}^{new} = argmax Q({\bf{\Theta, \Theta}}^{old})}$ 。
(4) 检查参数值的变化，如 ${||{\bf{\Theta}}^{old}-{\bf{\Theta}}^{new}||}$ ，若参数值收敛，算法结束；否则更新 ${{\bf{\Theta}}^{old}\gets{\bf{\Theta}}^{new}}$ ，返回第（2）步。

1.3.4 GMM的参数优化

理解了 EM 算法的基本原理后，我们现在就可以去推导 GMM 的参数优化流程了，从而证明式（1-9）和（1-10）的迭代方式是否有效。为了简便，这里只讨论一元的高斯分布，但同理也可以推广到多元的情况。

首先，一元 GMM 可表示为
1-26
根据式（1-25），我们构造需要优化的函数 ${Q({\bf{\Theta, \Theta}}^{old})}$ ，其中后验概率 ${p(y|x_i, {\bf{\Theta}}^{old})}$ 可根据式（1-9）计算，可得
1-27
前面我们一直假设 ${y}$ 的概率分布是已知的，但我们通常只能得到最终的输出 ${x}$ ，所以 ${\alpha_y}$ 实际上也是需要进行估计的参数，不过将其加到 ${\bf{\Theta}}$ 中同样是满足 EM 算法的要求的，为了方便这里还是对其单独讨论，因为从式（1-27）也可以看出 ${\alpha_y}$ 和 ${\bf{\Theta}}$ 极值点的计算是相互独立的。接下来，我们先对 ${y}$ 的概率分布进行估计。因为 ${y}$ 所有取值的概率之和为 1，所以这是一个等式约束的最优化问题。我们引入拉格朗日乘子 ${\lambda}$ ，那么其最优解满足
1-28
注意， ${p(y|x_i, {\bf{\Theta}}^{old})}$ 是基于旧参数计算的，所以在对新参数求导时相当于常数。那么求解式（1-28）可得
1-29
这样我们就可以知道，对于 ${Q({\bf{\Theta, \Theta}}^{old})}$ 来说， ${y}$ 等于某个值的概率的最优估计实际上就是各个样本在当前参数下属于该类的后验概率的平均值。同样，我们可以计算各个高斯分布的参数对于 ${Q({\bf{\Theta, \Theta}}^{old})}$ 的最优估计。因为
1-30
根据式（1-27）可得
1-31
同理可得
1-32
这样，我们就比较轻松地得到了对于 ${Q({\bf{\Theta, \Theta}}^{old})}$ 的各个高斯分布的参数的最优估计，而这刚好与我们最初猜想的式（1-10）是一致的，即首先根据现有参数计算出每个样本属于某个类的后验概率，然后以这个概率作为一种加权对该类的单高斯分布的参数进行估计，得到新的参数集合，并利用这个新的参数集进行迭代，从而在EM算法的理论支持下最终收敛到一个局部最优解，使得采样到的不完整数据的出现概率达到最大。

综上所述，使用 EM 算法求解 GMM 的参数的流程如图 1-4 所示。
图1-4

图1-4 EM 算法优化 GMM 参数流程

脚本示例

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

def get_gmm_data(n, means, covs, alpha):
    """ 
    生成 GMM 随机样本，这里采用了比较直观的方法，即首先根据 alpha 概率随机选择样本所属的高斯分布，
    然后再根据该高斯分布的参数生成随机生成最终的样本。这种方法相对比较低效，需要 n 次循环。

    Parameters
    ----------
    n : int
        生成随机样本的个数。
    means : list of 1D array
        各个高斯分布的均值，必须使用 list 来存放，或者使用 m x k 的二维矩阵，m 为高斯分布的个数。
    covs : list of 2D array
        各个高斯分布的协方差矩阵，必须使用 list 来存放，或者使用使用 m x k x k 的三维矩阵，k 为样本的维数。
    alpha : list or 1D array
        样本属于各个高斯分布的概率，总和必须为 1。
    
    Returns
    -------
    k x n 的二维数组，其中 k 为样本维数，可以从 means 推导。
    注意，为了方便进行矩阵的运算，脚本中所有样本都以列向量的形式存在，即每一列代表一个样本。
    """
    models = len(alpha)
    xs = []
    for i in range(n):
        j = np.random.choice(range(models), p=alpha)
        xs.append(np.random.multivariate_normal(means[j].reshape(-1), covs[j]))
    return np.array(xs).T

def gauss_prob(x, mean, cov):
    """
    根据给定 mean, cov 的高斯分布的概率密度函数计算样本 x 的概率密度。注意只适用于单高斯分布。

    Parameters
    ----------
    x : scalar or array
        可以为单个样本，也可以为多个样本。样本既可以是标量，也可以是列向量（单个样本也应该是 k x 1 的形式）。
    mean : scalar or array
        scalar 代表单变量的随机分布，array 代表多变量的随机分布。
        可以以行向量或列向量形式输入，list 形式也可以，实际上包括 scalar 都会首先被转换为 k x 1 的列向量形式。
    cov : scalar or 2D array
        如果 mean 是 scalar, cov 也必须为 scalar，否则应该为 k x k 的二维矩阵。
    
    Returns
    -------
    如果 x 是单个 scalar 样本，那么也返回一个 scalar，否则都返回一个 shape 为 (n,) 的 1D array。
    """
    return_scalar = True if np.shape(x) is () else False
    mean = np.reshape(mean, [-1, 1])
    dim = mean.shape[0]
    x = np.reshape(x, [dim, -1])
    cov = np.reshape(cov, [dim, dim])
    
    # 参考多维的高斯概率密度函数公式，注意这里可以一次性求得所有样本的概率密度
    tmp = (x - mean).T
    tmp = np.sum(tmp.dot(np.linalg.inv(cov)) * tmp, axis=1)
    prob = ((2*np.pi)**dim * np.linalg.det(cov))**(-0.5) * np.exp(-0.5 * tmp)
    
    if return_scalar:
        prob = prob[0]
    return prob

def gmm_em(xs, models):
    """
    基于 EM 算法估计 GMM 参数的具体实现。

    Parameters
    ----------
    xs : 2D array
        已有的样本，应该为 k x n 的形式，k 为样本维数，n 为样本数量。
        注意即便 k=1 样本也应该使用列向量的形式，因为代码大部分采用了矩阵运算。
    models : int
        高斯分布个数。
    
    Returns
    -------
    alpha : 1D array
        样本属于各个高斯分布的概率
    means : lsit of 1D array
        各个高斯分布的均值，以 list 形式存储
    covs : list of 2D array
        各个高斯分布的协方差矩阵，以 list 形式存储
    """
    # 根据样本的总体均值和协方差矩阵随机初始化各个高斯分布的参数，这种方式仅供参考
    dims, n = xs.shape[:2]
    mean = np.mean(xs, axis=1, keepdims=True)
    cov = (xs - mean).dot((xs - mean).T) / n
    alpha = np.random.rand(models) + 0.1
    alpha = alpha / np.sum(alpha)
    
    means, covs = [], []
    rnd_range = (np.max(xs, axis=1) - np.min(xs, axis=1)) * 0.8     # 样本数据的范围乘以一个缩放系数
    for i in range(models):
        means.append(mean.ravel() + (np.random.rand(dims) - 0.5) * rnd_range)   # 总体均值加上一个随机数
        covs.append(np.copy(cov))                                               # 直接使用总体协方差矩阵
    means, covs = np.array(means), np.array(covs)
    
    """ EM 算法流程 """
    for epoch in range(100):
        # 根据各个高斯分布的参数求得每个样本属于各个高斯分布的后验概率，即 E step
        probs = np.array([gauss_prob(xs, _mean, _cov) for _mean, _cov in zip(means, covs)])
        post_probs = np.diag(alpha).dot(probs)
        post_probs = post_probs / np.sum(post_probs, axis=0)

        # 根据求得的后验概率更新各个高斯分布的参数，即 M step
        new_alpha = np.sum(post_probs, axis=1)
        new_means = xs.dot(post_probs.T) / new_alpha
        new_covs = []
        for i in range(models):
            tmp = xs - new_means[:, [i]]
            new_covs.append(tmp.dot(np.diag(post_probs[i])).dot(tmp.T) / new_alpha[i])
        new_alpha = new_alpha / n
        new_means = new_means.T
        new_covs = np.array(new_covs)
        
        # 判断参数的更新是否停滞，这里使用了无穷阶范数，即 max(abs(x))，也可以改为 2nd 范数
        diff_alpha = np.linalg.norm(np.ravel(alpha - new_alpha), np.inf)
        diff_means = np.linalg.norm(np.ravel(means - new_means), np.inf)
        diff_covs = np.linalg.norm(np.ravel(covs - new_covs), np.inf)
        
        alpha = new_alpha
        means = new_means
        covs = new_covs
        
        if diff_alpha <= 1e-3 and diff_means <= 1e-3 and diff_covs <= 1e-3:
            print('Solved in {} epoches.'.format(epoch))
            break
    return alpha, list(means), list(covs)

if __name__ == '__main__':
    
    means = [np.array([5, 6]), np.array([1, 2])]
    covs = [np.array([[1, 0.1], [0.1, 0.5]]), np.array([[0.8, 1.2], [1.2, 3.5]])]
    alpha = np.array([0.3, 0.7])
    n = 2000
    
    xs = get_gmm_data(n, means, covs, alpha)
    #plt.scatter(xs[0], xs[1])
    #plt.show()
    
    _alpha, _means, _covs = gmm_em(xs, 2)
    print(_alpha)
    print(_means)
    print(_covs)

参考资料

http://www.seanborman.com/publications/EM_algorithm.pdf
https://statlect.com/fundamentals-of-probability/Jensen-inequality
https://statlect.com/fundamentals-of-probability/Kullback-Leibler-divergence