高斯混合模型_机器学习算法：高斯混合模型和EM算法

EM（Expectation-Maximization）算法是一种迭代式方法，主要应用于包含隐藏变量（latent variable）的参数估计，在无监督学习中有着广泛的应用，EM算法迭代包含两步：

• 利用估计的参数值来求对数似然期望（expectation）。
• 通过最大化对数似然期望来更新参数。

上述是EM算法的两个基本迭代部分，在实际应用中EM算法更多的看作是一种算法思想，而不是特定的算法步骤，所以接下来我将通过一些具体应用来进一步阐述EM算法的主要思想。

1、K-Means

K-Means算法是一种简单好用的聚类方法，之所以把它提到最前面来是因为我觉得同样作为无监督学习方法，K-Means在很多方面都深刻的体现着EM算法的思想。

那么首先我们来假设一个问题的形式：

input: 一组不带label的数据：

，现在要对其进行聚类。

output: 聚类模型

因为与监督学习不同，数据没有带label，所以对于之前介绍的逻辑回归等分类方法就无法应用，而是需要我们从数据自身进行发掘。下面给出K-Means算法的主要步骤：

1、随机初始化聚类重心（cluster centroids）:

2、Repeat until convergence：{
for every i, set:

for every j , set:

}

可以看到，K-Means主要包含两个迭代部分，首先是最小化重心与数据距离来更新参数，接着再对属于同一类的数据求平均值来更新重心。当进行了一段时间后重心不发生变化，K-Means算法也达到了收敛。

对应的，第一步更新参数类似于EM算法中的maximization，而更新中心则对应于expectation步。虽然和真正的EM算法有着很大的差别，但是这种迭代更新的思想却是不谋而合，所以我们可以将其看作一种简化版的EM来辅助理解。

2、高斯混合模型（GMM）

接下来要介绍的高斯混合模型（Gaussian Mixture Model）是EM算法的一个重要应用，它应用广泛，而很多时候EM算法是对GMM进行参数估计的有效手段。

2.1单高斯分布

高斯分布又名正态分布，可以说是我们最常见的一种分布之一了，而且实验表明许多受多个独立因素影响的数据（如体重身高等）都满足高斯分布。

如上图所示就是一个典型的二维高斯分布

高斯分布有两个重要参数：均值

以及方差

，若随机变量

服从高斯分布，则将其写作：

。

一维高斯分布的概率密度函数为：

多维高斯分布的概率密度函数为：

其中

是协方差矩阵的行列式。

2.2高斯混合模型

顾名思义，高斯混合模型就是多个高斯分布线性组合在一起形成的新的分布，理论上来说，任何分布都可以由多个混合分布组合表示，但由于高斯混合模型地应用广泛，这里我们主要讨论它。

如下图所示就是两个二维高斯分布的混合：

那么如何来描述上述的分布呢？很显然使用单个高斯模型是不够准确的（从图中可以直观地看到两个不同地聚类），所以由此引入了高斯混合模型地数学表达：

其中

是权重，且

，而

代表第k个分量模型。

如何理解GMM呢，首先我们先引入一个K维的二值变量

，其中某个特定的元素

，其余元素都为0，因此

满足

且

。

根据上述定义，

有K个不同的状态，用来表示观测数据

属于哪个分量模型，因此

每一个观测数据

都对应一个

，譬如

表示

来自第k个分量模型，同时变量

的边缘概率

则表示观测数据

属于第k个分量模型的概率。而由于

是“1 of K”的形式，因此变量

的分布也可以写成如下形似：

相应的，当对

指定一个特定的值时，那么

的条件概率分布就是一个特定的高斯分布：

也可以写成：

那么就可以得到联合概率分布

，而变量

的边缘概率分布就可以通过联合概率对所有可能的

值进行求和得到，即：

关于公式（2）第二个等式，乍一看会有些跳跃，但仔细想一想就会知道，所有可能的

值其实就是考虑

属于所有K个分量模型的可能性，所以体现的形式就是对所有模型加权求和。同时我们也发现经过一番推导我们得到了最初我们给出的高斯混合模型的表达式，所以我们所定义的

变量

可以看作是一个隐含变量

，用于决定观测数据属于哪个模型。

2.3使用EM方法估计GMM的参数

为了方便使用EM方法估计GMM的参数，我们首先定义一个新的量

，用它来表示后验概率

。根据贝叶斯理论，我们通过先验概率

以及条件概率

我们可以写出：

上面的概率公式给出的就是在已知观测样本的前提下估计隐藏变量的概率，这是我们在使用EM算法中的一个重要的过渡参数。

2.3.1极大似然函数

为了估计GMM的参数，我们要估计

这三个参数，所以我们首先将公式（1）的GMM模型改写成：

。然后由于

之间独立同分布，因此：

为了估计参数，我们要对其进行最大化似然函数操作，为了计算方便，我们首先给出公式（4）对数似然函数：

在最大化似然函数的前提下，我们首先来估计参数

，很自然的做法就是让公式（5）对变量

求导，然后令其导数为0，即：

接着对公式（6）两边同乘

后化简可得：

其中

，（化简借用公式（3））。

同样我们对

求导并令其导数为0，可得：

最后我们要估计的参数是

，但这里需要注意的是，在一开始我们对权值

有过限制，即：

，因此在极大化似然函数时要考虑这一条件。所以我们使用拉格朗日法乘数法来求下面式子的极大值：

然后可得：

对公式（10）前面分子分母约分后可得：

，同时又由于

，所以可得

。接着我们对公式（10）两边同乘

，进行化简后我们可以得到：

2.3.2 EM算法应用于GMM的算法步骤

根据（7）、（8）、（11）三式我们已经得到了GMM重要参数的更新迭代公式，这就是EM算法应用于高斯混合模型的实例，其过程包含E步（计算

）和M步（更新参数

）。接下来我们把具体的算法过程写下来：

Input：观测数据

：

（1）初始化GMM参数：均值

，协方差

，权值

（2） E步：根据当前模型的参数，计算 分模型k对观测数据

的响应度（responsibility）

（3） M步: 迭代更新模型的参数值：

其中：

（4）：估算对数似然函数的值：

重复（2）、（3）步直到算法收敛（似然函数值变化或参数变化小于设定值）

2.3.3 Example

下面我们使用一组真实数据来模拟使用EM算法求解GMM的过程。

数据：老忠实泉数据 （一组泉水喷发与时间的数据）

我们首先表示出数据的分布：

可以很明显地看出数据分为两类，因此我们可以使用包含两个高斯分布的GMM对其建模。

在实际操作时为了有利计算，往往会把数据进行归一化处理。

首先是对参数的初始化：

def init_params(shape, K): # K为模型中高斯分布数量
    N, D = shape # 样本数量和数据维度
    mu = np.random.rand(K, D)
    cov = np.array([np.eye(D)] * K)
    alpha = np.array([1.0 / K] * K)

    return mu, cov, alpha

然后是E步：

def getExpectation(Y, mu, cov, alpha):
    N = Y.shape[0]
    K = alpha.shape[0]

    # 计算响应度矩阵
    gamma = np.zeros((N, K))

    for k in range(K):
        gamma[:, k] = phi(Y, mu[k], cov[k])
    gamma = np.mat(gamma)

    # 计算每个模型对每个样本的响应度
    for k in range(K):
        gamma[:, k] = alpha[k] * gamma[:, k]
    for i in range(N):
        gamma[i, :] /= np.sum(gamma[i, :])
    return gamma

接着是M步：

def maximize(Y, gamma):
    N, D = Y.shape
    K = gamma.shape[1]

    #初始化参数值
    mu = np.zeros((K, D))
    cov = []
    alpha = np.zeros(K)

    # 更新每个模型的参数
    for k in range(K):
        # 第 k 个模型对所有样本的响应度之和
        Nk = np.sum(gamma[:, k])
        # 更新 mu
        # 对每个特征求均值
        for d in range(D):
            mu[k, d] = np.sum(np.multiply(gamma[:, k], Y[:, d])) / Nk
        print(mu)
        # 更新 cov
        cov_k = np.mat(np.zeros((D, D)))
        for i in range(N):
            cov_k += gamma[i, k] * (Y[i] - mu[k]).T * (Y[i] - mu[k]) / Nk
        cov.append(cov_k)
        # 更新 alpha
        alpha[k] = Nk / N
    cov = np.array(cov)
    return mu, cov, alpha

然后只需要重复迭代E步和M步达到收敛即可（可以自行设定收敛条件）：

def GMM(Y, K, times):
    Y = scale_data(Y)
    mu, cov, alpha = init_params(Y.shape, K)
    # 收敛条件
    delta = 0.000001
    for i in range(times):
        print("times: {}".format(i))
        gamma = Expectation(Y, mu, cov, alpha)
        last_mu = mu
        last_cov = cov
        last_alpha = alpha
        mu, cov, alpha = maximize(Y, gamma)
        mu_delta = np.sum(np.abs(last_mu - mu))
        cov_delta = np.sum(np.abs(last_cov - cov))
        alpha_delta = np.sum(np.abs(last_alpha - alpha))
        # 判断是否收敛
        if mu_delta < delta and cov_delta < delta and alpha_delta <delta:
            break

    return mu, cov, alpha

下面是迭代之后得到的结果：

从图中可以看到原始分布被有效的划分为了两个部分。

3、Conclusion

总结来说，EM算法的思路就是先求对数似然函数的期望值，然后再最大化似然函数以此来更新参数，通过不断地迭代来找的模型的参数。

对于使用EM算法比较重要的一点就是要找好隐含变量，通过隐含变量把观测数据与