聚類之高斯混合模型（Gaussian Mixture Model）

k-means應該是原來級別的聚類方法了，這整理下一個使用後驗機率準確評測其精度的方法—高斯混合模型。

咱們談到了用 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 的機率密度函數：

根據上面的式子，若是咱們要從 GMM 的分佈中隨機地取一個點的話，實際上能夠分爲兩步：首先隨機地在這 K 個 Component 之中選一個，每一個 Component 被選中的機率實際上就是它的係數 πk ，選中了 Component 以後，再單獨地考慮從這個 Component 的分佈中選取一個點就能夠了──這裏已經回到了普通的 Gaussian 分佈，轉化爲了已知的問題。

那麼如何用 GMM 來作 clustering 呢？其實很簡單，如今咱們有了數據，假定它們是由 GMM 生成出來的，那麼咱們只要根據數據推出 GMM 的機率分佈來就能夠了，而後 GMM 的 K 個 Component 實際上就對應了 K 個 cluster 了。根據數據來推算機率密度一般被稱做 density estimation ，特別地，當咱們在已知（或假定）了機率密度函數的形式，而要估計其中的參數的過程被稱做「參數估計」。

如今假設咱們有 N 個數據點，並假設它們服從某個分佈（記做 p(x) ），如今要肯定裏面的一些參數的值，例如，在 GMM 中，咱們就須要肯定 πk 、 μk 和 Σk 這些參數。 咱們的想法是，找到這樣一組參數，它所肯定的機率分佈生成這些給定的數據點的機率最大，而這個機率實際上就等於 ∏Ni=1p(xi) ，咱們把這個乘積稱做似然函數 (Likelihood Function)。一般單個點的機率都很小，許多很小的數字相乘起來在計算機裏很容易形成浮點數下溢，所以咱們一般會對其取對數，把乘積變成加和 ∑Ni=1logp(xi) ，獲得 log-likelihood function 。接下來咱們只要將這個函數最大化（一般的作法是求導並令導數等於零，而後解方程），亦即找到這樣一組參數值，它讓似然函數取得最大值，咱們就認爲這是最合適的參數，這樣就完成了參數估計的過程。

下面讓咱們來看一看 GMM 的 log-likelihood function ：

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

因爲在對數函數裏面又有加和，咱們無法直接用求導解方程的辦法直接求得最大值。爲了解決這個問題，咱們採起以前從 GMM 中隨機選點的辦法：分紅兩步，實際上也就相似於 K-means 的兩步。

• 估計數據由每一個 Component 生成的機率（並非每一個 Component 被選中的機率）：對於每一個數據 x_i 來講，它由第 k 個 Component 生成的機率爲
  
  γ(i,k)=πkN(xi|μk,Σk)∑Kj=1πjN(xi|μj,Σj)
  
  因爲式子裏的 μk和Σk 也是須要咱們估計的值，咱們採用迭代法，在計算 γ(i,k) 的時候咱們假定 μk和Σk 均已知，咱們將取上一次迭代所得的值（或者初始值）。
• 估計每一個 Component 的參數：如今咱們假設上一步中獲得的 γ(i,k) 就是正確的「數據 xi 由 Component k 生成的機率」，亦能夠當作該 Component 在生成這個數據上所作的貢獻，或者說，咱們能夠看做 xi 這個值其中有 γ(i,k)xi 這部分是由 Component k 所生成的。集中考慮全部的數據點，如今實際上能夠看做 Component 生成了 γ(1,k)x1,…,γ(N,k)xN 這些點。因爲每一個 Component 都是一個標準的 Gaussian 分佈，能夠很容易分佈求出最大似然所對應的參數值：
  
  μkΣk=1Nk∑i=1Nγ(i,k)xi=1Nk∑i=1Nγ(i,k)(xi−μk)(xi−μk)T
  
  其中 Nk=∑Ni=1γ(i,k) ，而且 πk 也瓜熟蒂落地能夠估計爲 Nk/N 。

重複迭代前面兩步，直到似然函數的值收斂爲止。
固然，上面給出的只是比較「直觀」的解釋，想看嚴格的推到過程的話，能夠參考 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×K 的矩陣，對於每個 xi ，咱們只要取該矩陣第 i 行中最大的那個機率值所對應的那個 Component 爲 xi 所屬的 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 自己只是一個模型，咱們這裏給出的迭代的辦法並非惟一的求解方法。感興趣的同窗能夠自行查找相關資料。

參考文獻：

• 大牛博客