漫谈 Clustering (3): Gaussian Mixture Model（转）

时间 2019-11-30

标签漫谈 clustering gaussian mixture model 繁體版

原文原文链接

上一次咱们谈到了用 k-means 进行聚类的方法，此次咱们来讲一下另外一个很流行的算法：Gaussian Mixture Model (GMM)。事实上，GMM 和 k-means 很像，不过 GMM 是学习出一些几率密度函数来（因此 GMM 除了用在 clustering 上以外，还常常被用于 density estimation ），简单地说，k-means 的结果是每一个数据点被 assign 到其中某一个 cluster 了，而 GMM 则给出这些数据点被 assign 到每一个 cluster 的几率，又称做 soft assignment 。算法

得出一个几率有不少好处，由于它的信息量比简单的一个结果要多，好比，我能够把这个几率转换为一个 score ，表示算法对本身得出的这个结果的把握。也许我能够对同一个任务，用多个方法获得结果，最后选取“把握”最大的那个结果；另外一个很常见的方法是在诸如疾病诊断之类的场所，机器对于那些很容易分辨的状况（患病或者不患病的几率很高）能够自动区分，而对于那种很难分辨的状况，好比，49% 的几率患病，51% 的几率正常，若是仅仅简单地使用 50% 的阈值将患者诊断为“正常”的话，风险是很是大的，所以，在机器对本身的结果把握很小的状况下，会“拒绝发表评论”，而把这个任务留给有经验的医生去解决。dom

废话说了一堆，不过，在回到 GMM 以前，咱们再稍微扯几句。咱们知道，无论是机器仍是人，学习的过程均可以看做是一种“概括”的过程，在概括的时候你须要有一些假设的前提条件，例如，当你被告知水里游的那个家伙是鱼以后，你使用“在一样的地方生活的是同一种东西”这相似的假设，概括出“在水里游的都是鱼”这样一个结论。固然这个过程是彻底“本能”的，若是不仔细去想，你也不会了解本身是怎样“认识鱼”的。另外一个值得注意的地方是这样的假设并不老是彻底正确的，甚至能够说老是会有这样那样的缺陷的，所以你有可能会把虾、龟、甚至是潜水员当作鱼。也许你以为能够经过修改前提假设来解决这个问题，例如，基于“生活在一样的地方而且穿着一样衣服的是同一种东西”这个假设，你得出结论：在水里有而且身上长有鳞片的是鱼。但是这样仍是有问题，由于有些没有长鳞片的鱼如今又被你排除在外了。机器学习

在这个问题上，机器学习面临着和人同样的问题，在机器学习中，一个学习算法也会有一个前提假设，这里被称做“概括偏执 (bias)”（bias 这个英文词在机器学习和统计里还有其余许多的意思）。例如线性回归，目的是要找一个函数尽量好地拟合给定的数据点，它的概括偏执就是“知足要求的函数必须是线性函数”。一个没有概括偏执的学习算法从某种意义上来讲毫无用处，就像一个彻底没有概括能力的人同样，在第一次看到鱼的时候有人告诉他那是鱼，下次看到另外一条鱼了，他并不知道那也是鱼，由于两条鱼总有一些地方不同的，或者就算是同一条鱼，在河里不一样的地方看到，或者只是看到的时间不同，也会被他认为是不一样的，由于他没法概括，没法提取主要矛盾、乎略次要因素，只好要求全部的条件都彻底同样──然而哲学家已经告诉过咱们了：世界上不会有任何样东西是彻底同样的，因此这我的即便是有无比强悍的记忆力，也绝学不到任何一点知识。函数

这个问题在机器学习中称做“过拟合 (Overfitting)”，例如前面的回归的问题，若是去掉“线性函数”这个概括偏执，由于对于 N 个点，咱们老是能够构造一个 N-1 次多项式函数，让它完美地穿过全部的这 N 个点，或者若是我用任何大于 N-1 次的多项式函数的话，我甚至能够构造出无穷多个知足条件的函数出来。若是假定特定领域里的问题所给定的数据个数老是有个上限的话，我能够取一个足够大的 N ，从而获得一个（或者无穷多个）“超级函数”，可以 fit 这个领域内全部的问题。然而这个（或者这无穷多个）“超级函数”有用吗？只要咱们注意到学习的目的（一般）不是解释现有的事物，而是从中概括出知识，并能应用到新的事物上，结果就显而易见了。学习

没有概括偏执或者概括偏执太宽泛会致使 Overfitting ，然而另外一个极端──限制过大的概括偏执也是有问题的：若是数据自己并非线性的，强行用线性函数去作回归一般并不能获得好结果。难点正在于在这之间寻找一个平衡点。不过人在这里相对于（如今的）机器来讲有一个很大的优点：人一般不会孤立地用某一个独立的系统和模型去处理问题，一我的天天都会从各个来源获取大量的信息，而且经过各类手段进行整合处理，概括所得的全部知识最终得以统一地存储起来，并能有机地组合起来去解决特定的问题。这里的“有机”这个词颇有意思，搞理论的人总能提出各类各样的模型，而且这些模型都有严格的理论基础保证能达到指望的目的，然而绝大多数模型都会有那么一些“参数”（例如 K-means 中的 k ），一般没有理论来讲明参数取哪一个值更好，而模型实际的效果却一般和参数是否取到最优值有很大的关系，我以为，在这里“有机”不妨看做是全部模型的参数已经自动地取到了最优值。另外，虽然进展不大，可是人们也一直都指望在计算机领域也创建起一个统一的知识系统（例如语意网就是这样一个尝试）。spa

废话终于说完了，回到 GMM 。按照咱们前面的讨论，做为一个流行的算法，GMM 确定有它本身的一个至关体面的概括偏执了。其实它的假设很是简单，顾名思义，Gaussian Mixture Model ，就是假设数据服从 Mixture Gaussian Distribution ，换句话说，数据能够看做是从数个 Gaussian Distribution 中生成出来的。实际上，咱们在 K-means 和 K-medoids 两篇文章中用到的那个例子就是由三个 Gaussian 分布从随机选取出来的。实际上，从中心极限定理能够看出，Gaussian 分布（也叫作正态 (Normal) 分布）这个假设实际上是比较合理的，除此以外，Gaussian 分布在计算上也有一些很好的性质，因此，虽然咱们能够用不一样的分布来随意地构造 XX Mixture Model ，可是仍是 GMM 最为流行。另外，Mixture Model 自己其实也是能够变得任意复杂的，经过增长 Model 的个数，咱们能够任意地逼近任何连续的几率密分布。scala

每一个 GMM 由个 Gaussian 分布组成，每一个 Gaussian 称为一个“Component”，这些 Component 线性加成在一块儿就组成了 GMM 的几率密度函数：code

$\displaystyle \begin{aligned} p(x) & = \sum_{k=1}^K p(k)p(x|k) \\ & = \sum_{k=1}^K \pi_k \mathcal{N}(x|\mu_k, \Sigma_k) \end{aligned}$

根据上面的式子，若是咱们要从 GMM 的分布中随机地取一个点的话，实际上能够分为两步：首先随机地在这个 Component 之中选一个，每一个 Component 被选中的几率实际上就是它的系数 $\pi_k$ ，选中了 Component 以后，再单独地考虑从这个 Component 的分布中选取一个点就能够了──这里已经回到了普通的 Gaussian 分布，转化为了已知的问题。component

那么如何用 GMM 来作 clustering 呢？其实很简单，如今咱们有了数据，假定它们是由 GMM 生成出来的，那么咱们只要根据数据推出 GMM 的几率分布来就能够了，而后 GMM 的个 Component 实际上就对应了个 cluster 了。根据数据来推算几率密度一般被称做 density estimation ，特别地，当咱们在已知（或假定）了几率密度函数的形式，而要估计其中的参数的过程被称做“参数估计”。orm

如今假设咱们有个数据点，并假设它们服从某个分布（记做），如今要肯定里面的一些参数的值，例如，在 GMM 中，咱们就须要肯定 $\pi_k$ 、 $\mu_k$ 和 $\Sigma_k$ 这些参数。咱们的想法是，找到这样一组参数，它所肯定的几率分布生成这些给定的数据点的几率最大，而这个几率实际上就等于 $\prod_{i=1}^N p(x_i)$ ，咱们把这个乘积称做似然函数 (Likelihood Function)。一般单个点的几率都很小，许多很小的数字相乘起来在计算机里很容易形成浮点数下溢，所以咱们一般会对其取对数，把乘积变成加和 $\sum_{i=1}^N \log p(x_i)$ ，获得 log-likelihood function 。接下来咱们只要将这个函数最大化（一般的作法是求导并令导数等于零，而后解方程），亦即找到这样一组参数值，它让似然函数取得最大值，咱们就认为这是最合适的参数，这样就完成了参数估计的过程。

下面让咱们来看一看 GMM 的 log-likelihood function ：

$\displaystyle \sum_{i=1}^N \log \left\{\sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)\right\}$

因为在对数函数里面又有加和，咱们无法直接用求导解方程的办法直接求得最大值。为了解决这个问题，咱们采起以前从 GMM 中随机选点的办法：分红两步，实际上也就相似于 K-means 的两步。

估计数据由每一个 Component 生成的几率（并非每一个 Component 被选中的几率）：对于每一个数据来讲，它由第个 Component 生成的几率为 $\displaystyle \gamma(i, k) = \frac{\pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)}{\sum_{j=1}^K \pi_j\mathcal{N}(x_i|\mu_j, \Sigma_j)}$
因为式子里的 $\mu_k$ 和 $\Sigma_k$ 也是须要咱们估计的值，咱们采用迭代法，在计算 $\gamma(i, k)$ 的时候咱们假定 $\mu_k$ 和 $\Sigma_k$ 均已知，咱们将取上一次迭代所得的值（或者初始值）。
估计每一个 Component 的参数：如今咱们假设上一步中获得的 $\gamma(i, k)$ 就是正确的“数据由 Component 生成的几率”，亦能够当作该 Component 在生成这个数据上所作的贡献，或者说，咱们能够看做这个值其中有 $\gamma(i, k)x_i$ 这部分是由 Component 所生成的。集中考虑全部的数据点，如今实际上能够看做 Component 生成了 $\gamma(1, k)x_1, \ldots, \gamma(N, k)x_N$ 这些点。因为每一个 Component 都是一个标准的 Gaussian 分布，能够很容易分布求出最大似然所对应的参数值： $\displaystyle \begin{aligned} \mu_k & = \frac{1}{N_k}\sum_{i=1}^N\gamma(i, k)x_i \\ \Sigma_k & = \frac{1}{N_k}\sum_{i=1}^N\gamma(i, k)(x_i-\mu_k)(x_i-\mu_k)^T \end{aligned}$
其中 $N_k = \sum_{i=1}^N \gamma(i, k)$ ，而且 $\pi_k$ 也瓜熟蒂落地能够估计为。
重复迭代前面两步，直到似然函数的值收敛为止。

固然，上面给出的只是比较“直观”的解释，想看严格的推到过程的话，能够参考 Pattern Recognition and Machine Learning 这本书的第九章。有了实际的步骤，再实现起来就很简单了。Matlab 代码以下：

（Update 2012.07.03：若是你直接把下面的代码拿去运行了，碰到 covariance 矩阵 singular 的状况，能够参见这篇文章。）

function varargout = gmm(X, K_or_centroids)
% ============================================================
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
%  - X: N-by-D data matrix.
%  - K_OR_CENTROIDS: either K indicating the number of
%       components or a K-by-D matrix indicating the
%       choosing of the initial K centroids.
%
%  - PX: N-by-K matrix indicating the probability of each
%       component generating each point.
%  - MODEL: a structure containing the parameters for a GMM:
%       MODEL.Miu: a K-by-D matrix.
%       MODEL.Sigma: a D-by-D-by-K matrix.
%       MODEL.Pi: a 1-by-K vector.
% ============================================================
 
    threshold = 1e-15;
    [N, D] = size(X);
 
    if isscalar(K_or_centroids)
        K = K_or_centroids;
        % randomly pick centroids
        rndp = randperm(N);
        centroids = X(rndp(1:K), :);
    else
        K = size(K_or_centroids, 1);
        centroids = K_or_centroids;
    end
 
    % initial values
    [pMiu pPi pSigma] = init_params();
 
    Lprev = -inf;
    while true
        Px = calc_prob();
 
        % new value for pGamma
        pGamma = Px .* repmat(pPi, N, 1);
        pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);
 
        % new value for parameters of each Component
        Nk = sum(pGamma, 1);
        pMiu = diag(1./Nk) * pGamma' * X;
        pPi = Nk/N;
        for kk = 1:K
            Xshift = X-repmat(pMiu(kk, :), N, 1);
            pSigma(:, :, kk) = (Xshift' * ...
                (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);
        end
 
        % check for convergence
        L = sum(log(Px*pPi'));
        if L-Lprev < threshold
            break;
        end
        Lprev = L;
    end
 
    if nargout == 1
        varargout = {Px};
    else
        model = [];
        model.Miu = pMiu;
        model.Sigma = pSigma;
        model.Pi = pPi;
        varargout = {Px, model};
    end
 
    function [pMiu pPi pSigma] = init_params()
        pMiu = centroids;
        pPi = zeros(1, K);
        pSigma = zeros(D, D, K);
 
        % hard assign x to each centroids
        distmat = repmat(sum(X.*X, 2), 1, K) + ...
            repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...
            2*X*pMiu';
        [dummy labels] = min(distmat, [], 2);
 
        for k=1:K
            Xk = X(labels == k, :);
            pPi(k) = size(Xk, 1)/N;
            pSigma(:, :, k) = cov(Xk);
        end
    end
 
    function Px = calc_prob()
        Px = zeros(N, K);
        for k = 1:K
            Xshift = X-repmat(pMiu(k, :), N, 1);
            inv_pSigma = inv(pSigma(:, :, k));
            tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);
            coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));
            Px(:, k) = coef * exp(-0.5*tmp);
        end
    end
end

函数返回的 Px 是一个 $N\times K$ 的矩阵，对于每个，咱们只要取该矩阵第行中最大的那个几率值所对应的那个 Component 为所属的 cluster 就能够实现一个完整的聚类方法了。对于最开始的那个例子，GMM 给出的结果以下：

相对于以前 K-means 给出的结果，这里的结果更好一些，左下角的比较稀疏的那个 cluster 有一些点跑得比较远了。固然，由于这个问题本来就是彻底有 Mixture Gaussian Distribution 生成的数据，GMM （若是能求得全局最优解的话）显然是能够对这个问题作到的最好的建模。

另外，从上面的分析中咱们能够看到 GMM 和 K-means 的迭代求解法其实很是类似（均可以追溯到 EM 算法，下一次会详细介绍），所以也有和 K-means 一样的问题──并不能保证老是能取到全局最优，若是运气比较差，取到很差的初始值，就有可能获得不好的结果。对于 K-means 的状况，咱们一般是重复必定次数而后取最好的结果，不过 GMM 每一次迭代的计算量比 K-means 要大许多，一个更流行的作法是先用 K-means （已经重复并取最优值了）获得一个粗略的结果，而后将其做为初值（只要将 K-means 所得的 centroids 传入 gmm 函数便可），再用 GMM 进行细致迭代。

如咱们最开始所讨论的，GMM 所得的结果（Px）不只仅是数据点的 label ，而包含了数据点标记为每一个 label 的几率，不少时候这其实是很是有用的信息。最后，须要指出的是，GMM 自己只是一个模型，咱们这里给出的迭代的办法并非惟一的求解方法。感兴趣的同窗能够自行查找相关资料。