码迷,mamicode.com
首页 > 其他好文 > 详细

GMM与K-means的那些事

时间:2015-07-25 18:31:26      阅读:444      评论:0      收藏:0      [点我收藏+]

标签:gmm   k-means   em   机器学习   


【原创】Liu_LongPo
转载请注明出处【CSDN】http://blog.csdn.net/llp1992

GMM算法

GMM ,Gaussian Mixture Model,顾名思义,就是说该算法由 多个高斯模型线性叠加 混合而成。每个高斯模型我们称之为 component 。GMM算法描述的是数据的本身存在的一种分布,如果component足够多的话,GMM可以逼近任意一种概率密度分布,但是此时的GMM模型就太过复杂了。

GMM 与 K-means

GMM算法跟K-means算法很像,它也是被用来聚类。在GMM中,component的个数我们定义为 K ,这跟K-means中的K一样,需要另外的方式去确定K值得大小。

在K-means中,每个聚类中心的初始化都会影响聚类效果,同样的,GMM算法对每个component的质心的初始化也很敏感。

解决方案

在K-means中,二分K-means算法就是为了弱化初始点对聚类效果的影响而提出来的,具体思路如下:

首先将所有点作为一个簇,然后将该簇一分为二。之后选择能最大程度降低聚类代价函数(也就是误差平方和)的簇划分为两个簇。以此进行下去,直到簇的数目等于用户给定的数目k为止

伪代码如下:

将所有数据点看成一个簇
当簇数目小于k时
对每一个簇
计算总误差
在给定的簇上面进行k-均值聚类(k=2)
计算将该簇一分为二后的总误差
选择使得误差最小的那个簇进行划分操作

同样的,GMM也可以有类似的二分GMM算法,关于二分K-means算法的具体流程以及代码请参考博客:

机器学习实战ByMatlab(四)二分K-means算法

GMM 模型

GMM由K个Gaussian分布线性叠加而成,先看看GMM的概率密度函数:

p(x)=k=1Kp(k)p(x|k)=k=1KπkN(x|μk,Σk)

该函数可以这么理解,假设我们有一个数据集,然后我们现在用GMM模型来描述这个数据集的分布。在已知数据集由component k 描述的情况下,数据集的概率密度函数为: p(x|k)

当然,总共有 K 个component,每个component 对生成数据集的贡献为 p(k) ,或者说数据集由component k生成的概率为 p(k),由 component k 生成数据集,其概率密度函数为 p(k)p(x|k) 。将所有的component加起来就得到了GMM的概率密度函数。有点绕口,大致懂就好。

如果我们要从 GMM 的分布中随机地取一个点的话,实际上可以分为两步:首先随机地在这 K个Gaussian Component 之中选一个,每个 Component 被选中的概率实际上就是它的系数 pi(k) ,选中了 Component 之后,再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布,转化为了已知的问题。

GMM 聚类

假设现在有一个大数据集,为什么要大数据集?待会会说。只要我们能用GMM算法来描述这个“客观存在”的数据集,那么GMM的K个component也就是对应的K个cluster了。根据数据来推算概率密度通常被称作 density estimation ,特别地,当我们在已知(或假定)了概率密度函数的形式,而要估计其中的参数的过程被称作“参数估计”。

每个component k 都是一个Gaussian分布,其均值设定为 μk ,方差设定为 Σk ,这个component的影响因子设定为 πk 。但是 我们一开始并不知道每个component k 的这几个参数的具体值,聚类误差函数中除了聚类后的label y之外,还有μkΣkπk这3个我们不知道的隐含变量,这时问题就得用EM算法来迭代求解。

参数与似然函数

现在假设我们有 N 个数据点,并假设它们服从某个分布(记作 p(x) ),现在要确定里面的一些参数的值,例如,在 GMM 中,我们就需要确定 影响因子 π(k)、各类均值 μk 和 各类协方差 Σk 这些参数。

我们的想法是,找到这样一组参数,它所确定的概率分布生成这些给定的数据点的概率最大,而这个概率实际上就等于 ΠNi=1p(xi),我们把这个乘积称作似然函数 (Likelihood Function)。通常单个点的概率都很小,许多很小的数字相乘起来在计算机里很容易造成浮点数下溢,因此我们通常会对其取对数,把乘积变成加和 Ni=1logp(xi) ,具体函数如下:

i=1Nlog{k=1KπkN(xi|μk,Σk)}

接下来我们只要将这个函数最大化(通常的做法是求导并令导数等于零,然后解方程),亦即找到这样一组参数值,它让似然函数取得最大值,我们就认为这是最合适的参数,这样就完成了参数估计的过程。

由于在对数函数里面又有加和,我们没法直接用求导解方程的办法直接求得最大值。为了解决这个问题,我们采取之前从 GMM 中随机选点的办法:分成两步,E步和M步,这就是用EM算法求解GMM的过程。其实这跟K-means的求解思想很像,或者说,K-means算法的求解中就是EM算法的精髓。

算法流程

1.估计数据由每个 Component 生成的概率(并不是每个 Component 被选中的概率):对于每个数据 来说,它由第 个 Component 生成的概率为:

γ(i,k)=πkN(xi|μk,Σk)Kj=1πjN(xi|μj,Σj)

其中 N(xi|μk,Σk) 就是后验概率

N(x|μ,Σ)=1(2π)D/21|Σ|1/2exp{?12(x?μ)TΣ?1(x?μ)}

2.通过极大似然估计可以通过求到令参数=0得到参数 μkΣk的值

μk=1Nki=1Nγ(i,k)xi

Σ=1Nki=1Nγ(i,k)(xi?μk)(xi?μk)T

其中,Nk=Ni=1γ(i,k) ,故 πk 可估计为 Nk/N

3.重复迭代前面两步,直到似然函数的值收敛为止。

Matlab 实现

GMM函数,来自 pluskid

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);
% isscalar 判断是否为标量
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
    %% Estiamtion Step
    Px = calc_prob();

    % new value for pGamma
    pGamma = Px .* repmat(pPi, N, 1);
    pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);

    %% Maximization Step
    % 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;  % 均值,也就是K类的中心
    pPi = zeros(1, K); % 概率
    pSigma = zeros(D, D, K); %协方差矩阵,每个都是 D*D

    % hard assign x to each centroids 
    % (X - pMiu)^2 = X^2 + pMiu^2 - 2*X*pMiu
    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)+diag(repmat(threshold,1,size(pSigma(:, :, k),1)))); % 方差矩阵求逆
        tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);
        coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma)); % det 求方差矩阵的行列式  
        Px(:, k) = coef * exp(-0.5*tmp);
    end
end
end

注意这句代码

coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma)); % det 求方差矩阵的行列式  

乍看之下好像跟我们的后验概率公式很不符合啊!后验概率公式如下:

N(x|μ,Σ)=1(2π)D/21|Σ|1/2exp{?12(x?μ)TΣ?1(x?μ)}

代码中的实现是 乘以协方差矩阵的逆矩阵的行列式的开根号后的值,根据线性代数的知识可以得到:|A?1|=1|A| , 故代码没错。

测试代码

clear all
clc

data = load(‘testSet.txt‘);
[PX,Model] = GMM(data,4);
[~,index] = max(PX‘);
cent = Model.Miu;
figure
I = find(index == 1);
scatter(data(I,1),data(I,2))
hold on
scatter(cent(1,1),cent(1,2),150,‘filled‘)
hold on
I = find(index == 2);
scatter(data(I,1),data(I,2))
hold on
scatter(cent(2,1),cent(2,2),150,‘filled‘)
hold on
I = find(index == 3);
scatter(data(I,1),data(I,2))
hold on
scatter(cent(3,1),cent(3,2),150,‘filled‘)
hold on
I = find(index == 4);
scatter(data(I,1),data(I,2))
hold on
scatter(cent(4,1),cent(4,2),150,‘filled‘)

效果如下

采用同样的数据集: GMM聚类效果如下:


技术分享

K-means效果如下:


技术分享

GMM与K-means对比

  • 聚类效果差不多,初始值敏感的缺点也差不多,解决方案也差不多,二分法。
  • 不同点:GMM输出的是数据点属于每个每类的概率,我们用最大似然方法去确定分类。就严谨性来说,用概率进行描述数据点的分类,GMM显然要比K-mean好很多。

关于GMM算法中奇异矩阵的问题

这个问题应该说是GMM算法的一个瑕疵,如果运行上面的代码可以发现在求协方差矩阵的逆矩阵的时候,有时候会报错,说协方差矩阵是奇异矩阵(singular),此时求解矩阵的行列式就会报错,关于这个问题,pluskid在 一篇博客
中已经说的很清楚。这里总结一下解决方法:

  1. 假设有N个D维的样本数据,如果协方差矩阵是奇异矩阵,那么要检查一下N是否足够大,D是否可以减少几维
  2. 最简单的解决方法之一,在协方差矩阵上加上一个对角线矩阵 λI,其中 λ 要足够小 ,不过如果数据集太小,这样也只是减少发生协方差矩阵是奇异矩阵的概率而已。

参考文献:漫谈 Clustering (3): Gaussian Mixture Model

版权声明:本文为博主原创文章,未经博主允许不得转载。

GMM与K-means的那些事

标签:gmm   k-means   em   机器学习   

原文地址:http://blog.csdn.net/llp1992/article/details/47058109

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!