統計學習方法ｃ++實現之八 EM算法與高斯混合模型

EM算法與高斯混合模型

前言

EM算法是一種用於含有隱變量的機率模型參數的極大似然估計的迭代算法。若是給定的機率模型的變量都是可觀測變量，那麼給定觀測數據後，就能夠根據極大似然估計來求出模型的參數，好比咱們假設拋硬幣的正面朝上的機率爲p（至關於咱們假設了機率模型），而後根據n次拋硬幣的結果就能夠估計出p的值，這種機率模型沒有隱變量，而書中的三個硬幣的問題（先拋A而後根據A的結果決定繼續拋B仍是C），這種問題中A的結果就是隱變量，咱們只有最後一個硬幣的結果，其中的隱變量沒法觀測，因此這種沒法直接根據觀測數據估計機率模型的參數，這時就須要對隱變量進行估計，進而獲得機率模型的參數，這裏要注意，機率模型是已知的（已經假定好了），包括隱變量的模型也是假設好的，只是具體的參數未知，這時候就須要用EM算法求解未知參數，這裏我用EM算法估計了高斯混合模型的參數，並用高斯混合模型實現了聚類，代碼地址。

EM算法

EM算法中文名稱是指望極大算法，EM是expectation maximization的縮寫，從名字就能夠窺視算法的核心，求指望，求極大。求誰的指望？求似然函數對隱變量的指望，因此，首先必須肯定隱變量是什麼。其次，對誰求極大？固然是求出機率模型的參數使得上一步的指望最大。算法以下：

輸入：觀測變量數據Y，隱變量數據Z（這裏也是知道的？其實這裏個人理解是，這裏不是已知的，可是倒是能夠根據假設的隱變量的參數獲得的），聯合分佈\(P(Y,Z|\theta)\), 條件分佈\(P(Z|Y,\theta)\)

輸出：模型參數\(\theta\)

1. E步:記\(\theta^{(i)}\)爲第i次迭代參數\(\theta\)的估計值，i+1次迭代的E步， 計算

   \(Q(\theta, \theta^{(i)})=E_Z[logP(Y,Z|\theta)|Y, \theta^{(i)}]=\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^{(i)})\)

2. M步，求使\(Q(\theta, \theta^{(i)})\)極大的\(\theta\),做爲下一次迭代的\(\theta^{(i+1)}\)

3. 重複2,3直到收斂

能夠看出最重要的在於求\(Q(\theta, \theta^{(i)})\)，那麼爲何每一次迭代最大化\(\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^{(i)})\)就能使觀測數據的似然函數最大（這是咱們的最終目的）？這裏書上有證實，很詳細，就不贅述了。先看一下咱們要最大化的極大似然函數，而後這裏主要引用西瓜書中的解釋來從理解EM算法：

\(L(\theta)=logP(Y|\theta)=log(\sum_ZP(Y,Z|\theta))=log(\sum_ZP(Y|Z,\theta)P(Z|\theta))\)

在迭代過程當中，若參數\(\theta\)已知，則能夠根據訓練數據推斷出最優隱變量Z的值（E步）；反之，若Z的值已知，則能夠方便的對參數\(\theta\)作極大似然估計(M步)。

能夠看出就是一種互相計算，一塊兒提高的過程。

這部分若是推導看不懂了，能夠結合下面的二維高斯混合模型的EM算法來理解。

c++實現

這裏我使用EM算法來估計高斯混合模型的參數來進行聚類，高斯混合模型還有一個很大的做用是進行前景提取，這裏僅僅用二維混合高斯模型進行聚類。

代碼結構

關鍵代碼

這裏實現起來其實沒什麼難點，難點在於推導參數的更新公式，詳情參考西瓜書p206。

void GMM::EMAlgorithm(vector<double> &alphaOld, vector<vector<vector<double>>> &sigmaOld,
        vector<vector<double>> &muOld) {
// compute gamma
    for (int i = 0; i < trainDataF.size(); ++i) {
        double probSum = 0;
        for (int l = 0; l < alpha.size(); ++l) {
            double gas = gaussian(muOld[l], sigmaOld[l], trainDataF[i]);
            probSum += alphaOld[l] * gas;
        }
        for (int k = 0; k < alpha.size(); ++k) {
            double gas = gaussian(muOld[k], sigmaOld[k], trainDataF[i]);
            gamma[i][k] = alphaOld[k] * gas / probSum;
        }
    }
// update mu, sigma, alpha
    for (int k = 0; k < alpha.size(); ++k) {
        vector<double> muNew;
        vector<vector<double>> sigmaNew;
        double alphaNew;
        vector<double> muNumerator;
        double sumGamma = 0.0;
        for (int i = 0; i < trainDataF.size(); ++i) {
            sumGamma += gamma[i][k];
            if (i==0) {
                muNumerator = gamma[i][k] * trainDataF[i];
            }
            else {
                muNumerator = muNumerator + gamma[i][k] * trainDataF[i];
            }
        }
        muNew = muNumerator / sumGamma;
        for (int i = 0; i < trainDataF.size(); ++i) {
            if (i==0) {
                auto temp1 = gamma[i][k]/ sumGamma * (trainDataF[i] - muNew);
                auto temp2 = trainDataF[i] - muNew;
                sigmaNew = vecMulVecToMat(temp1, temp2);
            }
            else {
                auto temp1 = gamma[i][k] / sumGamma * (trainDataF[i] - muNew);
                auto temp2 = trainDataF[i] - muNew;
                sigmaNew = sigmaNew + vecMulVecToMat(temp1, temp2);
            }
        }
        alphaNew = sumGamma / trainDataF.size();
        mu[k] = muNew;
        sigma[k] = sigmaNew;
        alpha[k] = alphaNew;
    }
}

總結

前面的代碼一直用vector來實現向量，可是這裏用到了矩陣，矩陣的相關計算都添加的計算函數。最正規的應該是寫個類，實現矩陣運算，可是這裏偷懶了，之後寫代碼必定要考慮周到，這樣添添補補的過低效了。