碼上快樂

相關內容    簡體    繁體

高斯混合聚類及EM實現

本文轉載自   查看原文   2016-06-04 17:49   2483    Machine Learning

一、引言

　　我們談到了用 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 由 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\}

　　在最大似然估計里面，由於我們的目的是把乘積的形式分解為求和的形式，即在等式的左右兩邊加上一個log函數，轉化為log后，還有log(a+b)的形式，因此，要進一步求解。

由於在對數函數里面又有加和，我們沒法直接用求導解方程的辦法直接求得最大值。為了解決這個問題，我們采取之前從 GMM 中隨機選點的辦法：分成兩步，實際上也就類似於 K-means 的兩步。

四、EM實現

  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. 重復迭代前面兩步，直到似然函數的值收斂為止。
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 就可以實現一個完整的聚類方法了。對於最開始的那個例子。

　　一個比較好的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 來最小化 J 並不容易，不過我們可以采取迭代的辦法：先固定 \mu_k ，選擇最優的 r_{nk}，很容易看出，只要將數據點歸類到離他最近的那個中心就能保證 J 最小。下一步則固定 r_{nk}，再求最優的 \mu_k。將 J 對 \mu_k 求導並令導數等於零，很容易得到 J 最小的時候 \mu_k 應該滿足：

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

亦即 \mu_k 的值應當是所有 cluster k 中的數據點的平均值。由於每一次迭代都是取到 J 的最小值，因此 J 只會不斷地減小（或者不變），而不會增加，這保證了 k-means 最終會到達一個極小值。雖然 k-means 並不能保證總是能得到全局最優解，但是對於這樣的問題，像 k-means 這種復雜度的算法，這樣的結果已經是很不錯的了。

下面我們來總結一下 k-means 算法的具體步驟：

  1. 選定 K 個中心 \mu_k 的初值。這個過程通常是針對具體的問題有一些啟發式的選取方法，或者大多數情況下采用隨機選取的辦法。因為前面說過 k-means 並不能保證全局最優，而是否能收斂到全局最優解其實和初值的選取有很大的關系，所以有時候我們會多次選取初值跑 k-means ，並取其中最好的一次結果。
  2. 將每個數據點歸類到離它最近的那個中心點所代表的 cluster 中。
  3. 用公式 \mu_k = \frac{1}{N_k}\sum_{j\in\text{cluster}_k}x_j 計算出每個 cluster 的新的中心點。
  4. 重復第二步，一直到迭代了最大的步數或者前后的 J 的值相差小於一個閾值為止。

　　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


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 EM及高斯混合模型 高斯混合和EM算法 4. EM算法-高斯混合模型GMM詳細代碼實現 混合高斯模型（Mixtures of Gaussians）和EM算法 EM算法和高斯混合模型GMM介紹 混合高斯分布與 EM 算法 EM算法原理以及高斯混合模型實踐 EM算法求高斯混合模型參數預計——Python實現 EM算法求高斯混合模型參數預計——Python實現 5. EM算法-高斯混合模型GMM+Lasso
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM