高斯混合聚类及EM实现

一、引言

　　我们谈到了用 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% 的阈值将患者诊断为“正常”的话，风险是非常大的，因此，在机器对自己的结果把握很小的情况下，会“拒绝发表评论”，而把这个任务留给有经验的医生去解决。

二、归纳法

　　一系列数据N要求对他拟合，如果不作要求的话，可以用一个N-1次多项式来完美拟合这N个点，比如拉格朗日插值，牛顿插值等，或者如果不限制次数的话可以找到无穷个完美函数，但是往往咱们要求指数型或者线性，就是需要对信息有机结合，GMM就是这样，要求用高斯模型来拟合，当然了可以构造任意的混合模型，不过根据中心极限定理等高斯模型比较合适。另外，Mixture Model 本身其实也是可以变得任意复杂的，通过增加 Model 的个数，我们可以任意地逼近任何连续的概率密分布。

三、原理分析

　　混合高斯模型基于多变量正态分布。类gmdistribution通过使用EM算法来拟合数据，它基于各观测量计算各成分密度的后验概率。

高斯混合模型常用于聚类，通过选择成分最大化后验概率来完成聚类。与k-means聚类相似，高斯混合模型也使用迭代算法计算，最终收敛到局部最优。高斯混合模型在各类尺寸不同、聚类间有相关关系的的时候可能比k-means聚类更合适。使用高斯混合模型的聚类属于软聚类方法（一个观测量按概率属于各个类，而不是完全属于某个类），各点的后验概率提示了各数据点属于各个类的可能性。

　　从几何上讲，单高斯分布模型在二维空间应该近似于椭圆，在三维空间上近似于椭球。遗憾的是在很多分类问题中，属于同一类别的样本点并不满足“椭圆”分布的特性。这就引入了高斯混合模型。

　　通观整个高斯模型，他主要是由方差和均值两个参数决定，,对均值和方差的学习,采取不同的学习机制,将直接影响到模型的稳定性、精确性和收敛性。

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

$\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 分布，转化为了已知的问题。

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

现在假设我们有个数据点，并假设它们服从某个分布（记作 p(x) ），现在要确定里面的一些参数的值，例如，在 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\}$

　　在最大似然估计里面，由于我们的目的是把乘积的形式分解为求和的形式，即在等式的左右两边加上一个log函数，转化为log后，还有log(a+b)的形式，因此，要进一步求解。

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

四、EM实现

假设模型参数已知，求概率。估计数据由每个 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$ 也顺理成章地可以估计为。
重复迭代前面两步，直到似然函数的值收敛为止。

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$ 的矩阵，对于每一个 x_i ，我们只要取该矩阵第行中最大的那个概率值所对应的那个 Component 为 x_i 所属的 cluster 就可以实现一个完整的聚类方法了。对于最开始的那个例子。

　　一个比较好的Matlab实现代码。

　　http://www.mathworks.com/matlabcentral/fileexchange/26184-em-algorithm-for-gaussian-mixture-model--em-gmm-

　　EM的推导参考

　　http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html

五、图形化显示

　　1. 为了展示高斯混合聚类的过程，先利用mvnrnd函数产生一些二变量高斯分布仿真数据：

mvnrnd(mu1,sigma1,200)

　　这句生成了200*2的矩阵。

mu1 = [1 2];
sigma1 = [3 .2; .2 2];
mu2 = [-1 -2];
sigma2 = [2 0; 0 1];
X = [mvnrnd(mu1,sigma1,200);mvnrnd(mu2,sigma2,100)];
scatter(X(:,1),X(:,2),10,'ko')

　　2. 然后拟合两成分高斯混合分布。这里，准确的成分数是已知的，而在处理实际数据时就需要通过测试、比较多个成分的拟合结果来决定这一数量。

　　options = statset('Display','final');

　　gm = gmdistribution.fit(X,2,'Options',options);

结果为

　　49 iterations, log-likelihood = -1207.91

　　3. 绘制估计的两成分混合分布的概率密度等高线如下。可以看到，两变量正态成分是相互重叠的，但是它们的峰值不同，这表明数据有理由分成两个聚类：

　　hold on

　　ezcontour(@(x,y)pdf(gm,[x y]),[-8 6],[-8 6]);

　　hold off

　　4. 使用cluster函数将拟合的混合分布数据分离成两类：

　　idx = cluster(gm,X);

　　cluster1 = (idx == 1);

　　cluster2 = (idx == 2);

　　scatter(X(cluster1,1),X(cluster1,2),10,'r+');

　　hold on

　　scatter(X(cluster2,1),X(cluster2,2),10,'bo');

　　hold off

　　legend('Cluster 1','Cluster 2','Location','NW')

每个聚类与混合分布中的一个成分有关。cluster基于估计的后验概率对各点做分类，属于哪个类的后验概率大则属于哪个类。posterior函数可以返回这个后验概率值。

例如，绘制各点属于第一成分的后验概率如下

P = posterior(gm,X);

scatter(X(cluster1,1),X(cluster1,2),10,P(cluster1,1),'+')

hold on

scatter(X(cluster2,1),X(cluster2,2),10,P(cluster2,1),'o')

hold off

legend('Cluster 1','Cluster 2','Location','NW')

clrmap = jet(80); colormap(clrmap(9:72,:))

ylabel(colorbar,'Component 1 Posterior Probability')

另一种分类的方法是使用前面计算的后验概率做软聚类。各点可以被赋予属于各类的一个得分，这个得分就是简单的后验概率，描述了各点与各聚类原型的相似程度。这样，各点可以按在聚类中的得分排序：

[~,order] = sort(P(:,1));

plot(1:size(X,1),P(order,1),'r-',1:size(X,1),P(order,2),'b-');

legend({'Cluster 1 Score' 'Cluster 2 Score'},'location','NW');

ylabel('Cluster Membership Score');

xlabel('Point Ranking');

尽管在数据散点图中很难看出分类的好坏，但得分图可以更明显的看出拟合的分布很好地完成了数据聚类，在得分0.5附近的点很少，即大多数点都易于分开。

使用高斯混合分布的软聚类与模糊k-means聚类方法相似，后者也是赋予各点相对各类的一个得分，但它进一步假设各聚类形状上近似球形，尺寸大小上也近似相等。这相当于所有成分协方差矩阵相同的高斯混合分布。而gmdistribution函数允许你设定各成分不同的协方差，默认情况下是为每个成分估计一个分离的无约束的协方差矩阵；而如果设定估计一个公共的对角协方差矩阵，则就与k-means相近了：

gm2 = gmdistribution.fit(X,2,'CovType','Diagonal',...

'SharedCov',true);

上面对于协方差的选项与模糊k-means聚类相似，但更灵活，允许不同的变量有不等的方差。

计算软聚类的得分可以不用先算硬聚类，可以直接使用posterior或cluster函数直接计算，如下

P2 = posterior(gm2,X); % equivalently [idx,P2] = cluster(gm2,X)

[~,order] = sort(P2(:,1));

plot(1:size(X,1),P2(order,1),'r-',1:size(X,1),P2(order,2),'b-');

legend({'Cluster 1 Score' 'Cluster 2 Score'},'location','NW');

ylabel('Cluster Membership Score');

xlabel('Point Ranking');

前面的例子中，混合分布的数据拟合与数据聚类是分开的，两步中使用了相同的数据。当然你也可以在cluster函数做聚类时使用新的数据点，实现新来点在原数据聚类中的分类。

　　1. 给定数据集X，首先拟合高斯混合分布，前面的代码已经实现

gm =

Gaussian mixture distribution with 2 components in 2 dimensions

Component 1:

Mixing proportion: 0.312592

Mean: -0.9082 -2.1109

Component 2:

Mixing proportion: 0.687408

Mean: 0.9532 1.8940

　　2. 然后可以对新数据集中的点分类到对X的聚类中

Y = [mvnrnd(mu1,sigma1,50);mvnrnd(mu2,sigma2,25)];

idx = cluster(gm,Y);

cluster1 = (idx == 1);

cluster2 = (idx == 2);

scatter(Y(cluster1,1),Y(cluster1,2),10,'r+');

hold on

scatter(Y(cluster2,1),Y(cluster2,2),10,'bo');

hold off

legend('Class 1','Class 2','Location','NW')

与前面的例子一样，各点的后验概率被作为得分，而不是直接做硬聚类。

完成聚类的新数据集Y应该来自与X相同的种群，即来自X的混合分布中，因为在估计Y的后验概率时使用了基于X拟合的混合高斯分布的联合概率。

六、讨论

　　6.1 kmeans

　　让我们先来看看它对于需要进行聚类的数据的一个基本假设吧：对于每一个 cluster ，我们可以选出一个中心点 (center) ，使得该 cluster 中的所有的点到该中心点的距离小于到其他 cluster 的中心的距离。虽然实际情况中得到的数据并不能保证总是满足这样的约束，但这通常已经是我们所能达到的最好的结果，而那些误差通常是固有存在的或者问题本身的不可分性造成的。例如下图所示的两个高斯分布，从两个分布中随机地抽取一些数据点出来，混杂到一起，现在要让你将这些混杂在一起的数据点按照它们被生成的那个分布分开来：

gaussian

　　由于这两个分布本身有很大一部分重叠在一起了，例如，对于数据点 2.5 来说，它由两个分布产生的概率都是相等的，你所做的只能是一个猜测；稍微好一点的情况是 2 ，通常我们会将它归类为左边的那个分布，因为概率大一些，然而此时它由右边的分布生成的概率仍然是比较大的，我们仍然有不小的几率会猜错。而整个阴影部分是我们所能达到的最小的猜错的概率，这来自于问题本身的不可分性，无法避免。因此，我们将 k-means 所依赖的这个假设看作是合理的。

　　基于这样一个假设，我们再来导出 k-means 所要优化的目标函数：设我们一共有 N 个数据点需要分为 K 个 cluster ，k-means 要做的就是最小化

$\displaystyle J = \sum_{n=1}^N\sum_{k=1}^K r_{nk} \|x_n-\mu_k\|^2$

　　这个函数，其中 $r_{nk}$ 在数据点 n 被归类到 cluster k 的时候为 1 ，否则为 0 。直接寻找 $r_{nk}$ 和 $\mu_k$ 来最小化并不容易，不过我们可以采取迭代的办法：先固定 $\mu_k$ ，选择最优的 $r_{nk}$ ，很容易看出，只要将数据点归类到离他最近的那个中心就能保证最小。下一步则固定 $r_{nk}$ ，再求最优的 $\mu_k$ 。将对 $\mu_k$ 求导并令导数等于零，很容易得到最小的时候 $\mu_k$ 应该满足：

$\displaystyle \mu_k=\frac{\sum_n r_{nk}x_n}{\sum_n r_{nk}}$

亦即 $\mu_k$ 的值应当是所有 cluster k 中的数据点的平均值。由于每一次迭代都是取到的最小值，因此只会不断地减小（或者不变），而不会增加，这保证了 k-means 最终会到达一个极小值。虽然 k-means 并不能保证总是能得到全局最优解，但是对于这样的问题，像 k-means 这种复杂度的算法，这样的结果已经是很不错的了。

下面我们来总结一下 k-means 算法的具体步骤：

选定 K 个中心 $\mu_k$ 的初值。这个过程通常是针对具体的问题有一些启发式的选取方法，或者大多数情况下采用随机选取的办法。因为前面说过 k-means 并不能保证全局最优，而是否能收敛到全局最优解其实和初值的选取有很大的关系，所以有时候我们会多次选取初值跑 k-means ，并取其中最好的一次结果。
将每个数据点归类到离它最近的那个中心点所代表的 cluster 中。
用公式 $\mu_k = \frac{1}{N_k}\sum_{j\in\text{cluster}_k}x_j$ 计算出每个 cluster 的新的中心点。
重复第二步，一直到迭代了最大的步数或者前后的的值相差小于一个阈值为止。

　　6.2 GMM与kmeans的区别

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

　　如我们最开始所讨论的，GMM 所得的结果（Px）不仅仅是数据点的 label ，而包含了数据点标记为每个 label 的概率，很多时候这实际上是非常有用的信息。最后，需要指出的是，GMM 本身只是一个模型，我们这里给出的迭代的办法并不是唯一的求解方法。感兴趣的同学可以自行查找相关资料。

GMM中几个高斯模型就是几个聚类簇，也就是说不能聚类是K个，我想用M个或者任意多个高斯模型，而且假设簇终点服从高斯分布，因此个人认为FCM使用更广泛。

http://blog.pluskid.org/?p=39

http://www.cnblogs.com/CBDoctor/archive/2011/11/06/2236286.html

http://www.cnblogs.com/zhangchaoyang/articles/2624882.html

http://chunqiu.blog.ustc.edu.cn/?p=437

posted @ 2016-06-04 17:49 加拿大小哥哥阅读(6847) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

高斯混合聚类及EM实现

一、引言

二、归纳法

三、原理分析

四、EM实现

五、图形化显示

六、讨论

6.1 kmeans

6.2 GMM与kmeans的区别

公告

　　6.1 kmeans

　　6.2 GMM与kmeans的区别