Categories

Calendar

January 2022
M T W T F S S
« Jun    
 12
3456789
10111213141516
17181920212223
24252627282930
31  

漫谈 Clustering (3): Gaussian Mixture Model

cluster_logo本文是“漫谈 Clustering 系列”中的第 4 篇,参见本系列的其他文章。

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

废话说了一堆,不过,在回到 GMM 之前,我们再稍微扯几句。我们知道,不管是机器还是人,学习的过程都可以看作是一种“归纳”的过程,在归纳的时候你需要有一些假设的前提条件,例如,当你被告知水里游的那个家伙是鱼之后,你使用“在同样的地方生活的是同一种东西”这类似的假设,归纳出“在水里游的都是鱼”这样一个结论。当然这个过程是完全“本能”的,如果不仔细去想,你也不会了解自己是怎样“认识鱼”的。另一个值得注意的地方是这样的假设并不总是完全正确的,甚至可以说总是会有这样那样的缺陷的,因此你有可能会把虾、龟、甚至是潜水员当做鱼。也许你觉得可以通过修改前提假设来解决这个问题,例如,基于“生活在同样的地方并且穿着同样衣服的是同一种东西”这个假设,你得出结论:在水里有并且身上长有鳞片的是鱼。可是这样还是有问题,因为有些没有长鳞片的鱼现在又被你排除在外了。

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

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

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

废话终于说完了,回到 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 的个数,我们可以任意地逼近任何连续的概率密分布。

每个 GMM 由 K 个 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 的分布中随机地取一个点的话,实际上可以分为两步:首先随机地在这 K 个 Component 之中选一个,每个 Component 被选中的概率实际上就是它的系数 \pi_k ,选中了 Component 之后,再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布,转化为了已知的问题。

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

现在假设我们有 N 个数据点,并假设它们服从某个分布(记作 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\}

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

  1. 估计数据由每个 Component 生成的概率(并不是每个 Component 被选中的概率):对于每个数据 x_i 来说,它由第 k 个 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 均已知,我们将取上一次迭代所得的值(或者初始值)。

  2. 估计每个 Component 的参数:现在我们假设上一步中得到的 \gamma(i, k) 就是正确的“数据 x_i 由 Component k 生成的概率”,亦可以当做该 Component 在生成这个数据上所做的贡献,或者说,我们可以看作 x_i 这个值其中有 \gamma(i, k)x_i 这部分是由 Component k 所生成的。集中考虑所有的数据点,现在实际上可以看作 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 也顺理成章地可以估计为 N_k/N 。

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

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

(Update 2012.07.03:如果你直接把下面的代码拿去运行了,碰到 covariance 矩阵 singular 的情况,可以参见这篇文章。)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
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 ,我们只要取该矩阵第 i 行中最大的那个概率值所对应的那个 Component 为 x_i 所属的 cluster 就可以实现一个完整的聚类方法了。对于最开始的那个例子,GMM 给出的结果如下:

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 本身只是一个模型,我们这里给出的迭代的办法并不是唯一的求解方法。感兴趣的同学可以自行查找相关资料。

199 comments to 漫谈 Clustering (3): Gaussian Mixture Model

  • bendy

    两个Gaussian分布相加应该是Gaussian分布。您说的“在平面上,要将两个高斯分布的数据类别分开来的最优的分界线是二次曲线,如果两个分布的 variance 是一样的话,还可以变成简单的直线。”这一性质,我到没见过,没有太大的概念。是不是应该在数值分析后可以显示出来?

    • 就是两个 Gaussian 相减等于零(就是分隔的那个曲线/面),这个解析式是一个二次的,印象中《Pattern Classification》好像有详细分析三维中的各种情况的。

  • summer_night

    soft assignment 和GMM 应该不是完全相等吧,他们有什么区别和联系呢?

  • Friday

    我想知道K值如何估算?还有πk 如何估算?我是初学者~请指点,谢谢

  • Friday

    师兄研二啊,今年?

  • Friday

    假设我有一批数据,希望用GMM进行拟合,那么我该如何操作?这个K的给定有什么方法或者技巧?非常感谢

    • 根据先验知识来定,或者枚举找最优的吧。也有方法可以自动估计类别数目的,比如 Data Spectroscopy: Learning Mixture Models using Eigenspaces of Convolution Operators 。不过没有尝试过,不知道在实际数据上效果如何。

  • Friday

    我又在文献中看到过hold-out方法,就是枚举吧?从K=1开始,一直到K=n+1时过拟合就停止,将K=n视为optimal,是不是可以?
    还有如果是用K-mean的结果作为初值,那么k-mean聚类出的cluster数目是不是可以作为K的初值?然后其他的参数一并计算?

  • […] overfit å’Œ underfit 的问题可以参见这篇文章的 3~6 段,这位学长讲得非常深入浅出了。课上提到的一种解决方法叫 […]

  • Friday

    lz,请问:为了决定最优的K值,我先用K-mean,从K=1开始进行聚类,然后K+1……依次进行下去,直到K+1时的J大于K时的J,即过拟合,就停止,选取此时的K作为GMM的K,可以么?

    • ä½ çš„ J 是 kmeans 的目标函数吗?如果是那样的话,是行不通的,因为 kmeans 的目标函数肯定随着 k 增大而减小,当 k 等于数据点个数的时候 J = 0 。如果想要自动确定 GMM 类别个数的话,确实是有一些工作的,比如 Belkin çš„ Data Spectroscopy: Learning Mixture Models using Eigenspaces of Convolution Operators 这个工作。实际我没有试过,不知道在实际数据中效果如何。

  • Friday

    那该怎么确定GMM的类别个数啊?
    还有一个疑问,这个GMM做出来的应该是连续的,而且它的每个部分也是连续的;那么你在找最优的参数的时候,怎么能计算一个点的概率?不是没有意义么?
    辛苦楼主~

  • Friday

    我按照你的代码执行了一遍,为什么这个矩阵Px里面有的值大小超过了1?这个是可以的??

    • 你好,Px 只是把数据点 x 带入到 Gaussian model 里得到的值,可以超过 1 的。

      • Friday

        lz,假设一下:一堆状态点,每个点是一个三维的数据,我用了GMM建模,找到了其联合pdf;现在我有一个新的状态点,我要求它的概率,显然单点的是没有意义的,而且在小区域内的概率也是很小的吧,那么这个除了将点归类还有什么其他的意义?这个问题对我来说很棘手,请楼主原谅问题的愚蠢~谢谢~~·

        • 你好,不好意思,我不太明白你想要解决的问题是什么,或者你想要做什么。

  • 姜

    你好,我现在在写ICP方面的论文,看到Gaussian Mixture Models后,感觉可以在这方面下手,提出一种新的方法,请问,您能给出解决散文点云空间分布的高斯混合模型的C++代码吗?非常感谢!!

  • Friday

    我是想求某个状态点的概率

    • 像 Gaussian 分布这种绝对连续分布单点的概率是为零的。

      • Friday

        我知道……您见我的邮件吧……我实在描述不清楚了,麻烦您了

  • ke

    非常棒的文章。

    另外BOA就相当于美国的建行 : D

    我们实验室也好几个浙大的师兄,都是非常棒,祝贺楼主拿到MIT的offer(我是不是回错地方了,不过这两天一直琢磨您这篇文章来的,快考试了= =|||)

  • eagler_hu

    能否推荐两本讲聚类算法或模式识别方面较经典的书籍?!谢谢!

  • 文章很强大

  • 清乐时啸

    受教
    转载一下,谢谢!
    如不宜,请email通知我。

  • xiaotao

    师兄好,能否给解释一下高斯模型公式中T是指的什么?就是exp的次幂里的那个T

  • xiaotao

    师“熊”的文章果然好强大!!看了确实受益匪浅……我依然有两个小问题不明:
    1)我学的是用混合高斯模型对背景建模方面的,假设一个像素点与模型第一次匹配成功后,就可以说这个像素点是背景点了,那么还要继续和后面的(2至4个)模型匹配吗?
    2)更新后权值怎么保证累加和为1呢?或者说权值累加和必须要等于1吗?
    3)alpha的值是自己设定的经验值吗?大概为多少?
    跪谢!!!

  • Friday

    师兄,请问在matlab运行到第96行( inv_pSigma = inv(pSigma(:, :, k));)时,警告说矩阵奇异,为什么会出现这种情况啊?按道理来说,对于聚类的数据应该是没有特殊要求的啊!
    非常感谢!!打扰了!

    • 你好,这个 code 主要是 demo code ,为了结构清晰很多地方没有做细致处理,例如那里求逆的时候我们一般是不会直接用 inv 来求的,可以用 \ 或者 / 直接来解方程。另外这里也没有对 covariance 做任何限制或者 regularization ,如果数据量不够多的话这是不好的并且很可能出现你说的这种情况。所以,如果你是想跑实验的话,最好是去找网上其他专门写来跑实验的 GMM 代码,不仅处理了各种情况,而且速度也会快许多。 :)

      • Friday

        麻烦师兄给个程序……我没找到啊……谢谢了

  • DJ

    哈哈,还是kid的blog好

  • GMM

    文章太强大了,看不太懂呢!师兄,恕我无知,你那个程序运行出来的结果varargout应该是什么样的矩阵?我想输入一幅高清图像的灰度值矩阵,为21025*200的矩阵,给其分15类,但貌似总是运行不了,学长能解释一下吗?还有,第83行 [dummy labels] = min(distmat, [], 2)中的dummy是做什么用的,貌似出错。。。

  • 背景GMM

    代码中将每个数据点看做是D维的,如果将NxD数据替换为一幅图像,可以直接用吗? % new value for pGamma
    pGamma = Px .* repmat(pPi, N, 1);
    pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);
    第一行中pGamma 不已经是概率了吗?为什么还要第二行?请教一下!先谢了

    • Florida

      正好读了下,楼主分析的很好很强大~顺便回答下这位童鞋

      第一行是求的step 1中公式的分子,还没有计算完;第二行是除以公式中的分母,结果才是 soft membership~ lol~

  • llxlf2012

    那么GMM是计算将数据点assign到某个类别的概率,这个作用和FCM模糊聚类有什么差异?

  • […] 我之前写过一篇介绍 Gaussian Mixture Model (GMM) 的文章,并在文章里贴了一段 GMM 实现的 Matlab 示例代码,然后就不断地有人来问我关于那段代码的问题,问得最多的就是大家经常发现在跑那段代码的时候估计出来的 Covariance Matrix 是 singular 的,所以在第 96 行求逆的时候会挂掉。这是今天要介绍的主要话题,我会讲得罗嗦一点,把关于那篇文章我被问到的一些其他问题也都提到一下,不过,在步入正题之前,不妨先来小小地闲扯一下。 […]

  • […] a comment » 我之前写过一篇介绍 Gaussian Mixture Model (GMM) 的文章,并在文章里贴了一段 GMM 实现的 Matlab […]

  • […] 如果说 K-means 和 GMM 这些聚类的方法是古代流行的算法的话,那么这次要讲的 Spectral Clustering 就可以算是现代流行的算法了,中文通常称为“谱聚类”。由于使用的矩阵的细微差别,谱聚类实际上可以说是一“类”算法。 […]

  • gs

    你好:请教一个问题,多维的数据,如何可视化聚类结果!就是上面那个图是怎么画的?

  • Everest

    第三次来看这篇文章,前两次看了不太明白,回去查了不少资料,今天终于全明白了…感谢楼主~!

  • Alston Tsao

    你好,小弟今天發現您的blog,拜讀了這一整個系列的文章,以及其他許多,受益良多。很少遇到像你這樣能夠把數學解釋的清楚透徹並且附上source code的,無論是python or matlab,把source code實際載下來玩能夠更深刻地體會抽象的演算法,非常感謝您無私的分享,學習了 lol。

  • 晶晶love喵咪

    师兄,您好,我是研一的新生,您的这篇文章我已经看了好久了,并且基于你的这篇文章的代码,重新编写了一些代码,也是关于高斯混合模型的,只不过是将每一个component看成是由一个高斯混合模型组成的,也就是将原先的高斯分布假设成了高斯混合分布,

  • 晶晶love喵咪

    然后就是在求高斯分布时的协方差矩阵求逆时会出现奇异的情况,不知道师兄可否给点建议,鄙人将不甚感激~

    • 你是说,两层的混合高斯模型?这个和普通的混合高斯模型不是一样的吗?

    • AdamD

      如果使用EM-MAP聚类的话,一个先验假设是每个component具有相同的权值,也就是Mixture里面那个alpha,这个假设是不是有点arbitrary,不知道之前有没有人做过优化这个权值的工作呢?

  • 晶晶love喵咪

    师兄说的两层高斯混合模型就是把类用高斯混合模型拟合吧,原先高斯模型是用高斯分布拟合的,应该就是我做的模型,我是基于半监督学习上的高斯混合模型分类,函数输进去的是已知数据的先验概率和各类高斯混合模型的中心点坐标。更新的值是先验概率和中心点坐标的值,采取迭代的方式,不知道是我代码出错还是有什么限制,分出的图像效果不好,望师兄指点。

    • 我不是很清楚你是在做什么,就是每个类别分别 train 出 GMM ,然后分类的时候看每个类别的似然?GMM çš„ covariance singular 的话,可以加 regularization,具体可以看我这篇 blog : http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/ 。分类效果不好的话,这个可能有很多因素,我也不知道了呀。

  • Steven_Yang

    师兄真是热心人……
    厉害!

  • forrest1991

    为什么第k个协方差矩阵是pSigma(:, :, k)而不是pSigma(k, :,:),这样做有什么好处?

  • forrest1991

    您好,用pinv代替inv是不是就能处理pSigma不可逆的情况?

    • 你好,应该是可以解决计算的问题,它忽略掉那些零的特征值进行求逆的话,应该也可以从某种角度来进行 regularization 的解释,不过 pinv 计算开销挺大的。

  • liubinzyx

    你好,看了你的文章收获很大啊,很牛!!!但是我还有个问题想问一下,高斯混合模型中的各个component之间是否相互独立呢?

  • liubinzyx

    不是独立正态分布相加还是正态分布吗,那这些component相加应该还是一个正态分布?不知道我这样理解可不可以

    • “两个正太分布随机变量的和”和“两个正太分布概率密度函数的和”是不同的东西。

      • liubinzyx

        哦,看样子,还有些概念没弄清楚,谢谢指导啊,我再琢磨一下

  • jetfish

    Excellent!博主,请多谢文章,受益匪浅!!

  • jetfish

    博主,该怎么将结果绘图出来了?

  • jetfish

    急求回复,博主,谢谢,谢谢!!

  • jetfish

    博主,那怎么讲一个1000个点的坐标值转为N*D的X矩阵呀?我现在的数据是2*1000的矩阵。

  • lz,如何求解一维混合高斯分布的极值(点)?谢谢!

  • systolic

    很好的东西,讲得很透彻。
    在code的注释中,
    MODEL.Miu: a K-by-D matrix

    如果用一个D-by-K matrix,本质上是一样的,但物理意义上是否和另外两个参数更一致些?

  • AdamD

    如果使用EM-MAP聚类的话,一个先验假设是每个component具有相同的权值,也就是Mixture里面那个alpha,这个假设是不是有点arbitrary,不知道之前有没有人做过优化这个权值的工作呢?

    或者也许是这个问题不太重要?

  • sbshiu

    你好
    程式碼中的
    coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));

    det(inv_pSigma)是不是要用pSigma而不是inv_pSigma?

    感謝

  • Moon

    求问函数最后返回的时候 varargout = {Px, model}为什么是Px而不是 pGamma?
    Px应该只是高斯函数的值,pGamma才是后验概率啊。。。

    • 输出参数的后验概率来做什么?求的是数据的概率呀。当然你想要的话也可以输出来没有什么关系……

      • Moon

        不好意思是我弄错了- -。。。应该是比较数据点x属于各个特定分布的概率来决定x属于哪个簇,而不是求x由该分布生成的概率,多谢~

      • Moon

        不好意思,可能是我才疏学浅,仔细想想还是无法理解。。。《Pattern Recognition and Machine Learning》的9.3.2讲与K-means的关系的时候,第一段说“makes a soft assignment based on the posterior probabilities”,这个soft assignment应该就是r(z_{n,k}),也就是数据点被 assign 到每个 cluster 的后验概率啊。。。而K-means中数据点的hard ssignment也是根据r(z_{n,k})的极限值得到的,为什么这里不是根据比较r(z_{n,k})的大小来确定数据x属于哪个簇呢?

        • 哦,我明白你的意思了。两者的区别是 P(x|k) å’Œ P(k|x),最后归为 hard decision 的时候可以看成分别是 ML å’Œ MAP 两种方法,各自也是有道理的,不过如果要和 PRML 上的一致的话,应该确实是用你说的这种 MAP 做 hard decision çš„~

        • Moon

          明白啦,感谢~

  • will

    tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);
    这里是不是漏掉一个转置符号,在前一个Xshift上?

  • huda

    文章不错,赞一个,但是程序中
    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, :);
    不懂,跪求解答!

  • may

    我想问一下,混合高斯模型与多个服从正态分布的随机向量线性组合之间有什么关系?

  • LEEE

    您好!现在正看GMM,您的文章让我受益匪浅,感谢您的分享,谢谢。

Leave a Reply

 

 

 

You can use these HTML tags

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>