EM算法求高斯混合模型參數預計——Python實現

EM算法一般表述：

     

       當有部分數據缺失或者無法觀察到時，EM算法提供了一個高效的迭代程序用來計算這些數據的最大似然預計。在每一步迭代分為兩個步驟：期望（Expectation）步驟和最大化（Maximization）步驟。因此稱為EM算法。

       如果所有數據Z是由可觀測到的樣本X={X1, X2。……, Xn}和不可觀測到的樣本Z={Z1, Z2，……, Zn}組成的。則Y = X∪Z。EM算法通過搜尋使所有數據的似然函數Log(L(Z; h))的期望值最大來尋找極大似然預計。注意此處的h不是一個變量，而是多個變量組成的參數集合。此期望值是在Z所遵循的概率分布上計算，此分布由未知參數h確定。然而Z所遵循的分布是未知的。EM算法使用其當前的如果h`取代實際參數h，以預計Z的分布。

                                                             Q( h`| h) = E [ ln P(Y|h`) | h, X ]

       EM算法反復下面兩個步驟直至收斂。

       步驟1：預計（E）步驟：使用當前如果h和觀察到的數據X來預計Y上的概率分布以計算Q( h` | h )。

                                                             Q( h` | h ) ←E[ ln P(Y|h`) | h, X ]

       步驟2：最大化（M）步驟：將如果h替換為使Q函數最大化的如果h`:

                                                              h ←argmaxQ( h` | h )

高斯混合模型參數預計問題：

          簡單起見，本問題研究兩個高斯混合模型參數預計k=2。

       問題描寫敘述：如果X是由k個高斯分布均勻混合而成的。這k個高斯分布的均值不同，可是具有同樣的方差。

設樣本值為x1, x2, ……, xn。xi能夠表示為一個K+1元組< xi, zi1, zi2, …, zik>，當中僅僅有一個取1，其余的為0。

此處的zi1到zik為隱藏變量。是未知的。且隨意zij被選擇的概率相等。即

                                                 P（zij = 1）=1/k (j=1,2,3.....k)

       EM算法求解過程推導例如以下：

   

Python實現（模擬2個正態分布的均值預計）：

#coding:gbk
import math
import copy
import numpy as np
import matplotlib.pyplot as plt

isdebug = False

# 指定k個高斯分布參數。這里指定k=2。注意2個高斯分布具有同樣均方差Sigma，分別為Mu1,Mu2。

def ini_data(Sigma,Mu1,Mu2,k,N):
    global X
    global Mu
    global Expectations
    X = np.zeros((1,N))
    Mu = np.random.random(2)
    Expectations = np.zeros((N,k))
    for i in xrange(0,N):
        if np.random.random(1) > 0.5:
            X[0,i] = np.random.normal()*Sigma + Mu1
        else:
            X[0,i] = np.random.normal()*Sigma + Mu2
    if isdebug:
        print "***********"
        print u"初始觀測數據X："
        print X
# EM算法：步驟1，計算E[zij]
def e_step(Sigma,k,N):
    global Expectations
    global Mu
    global X
    for i in xrange(0,N):
        Denom = 0
        for j in xrange(0,k):
            Denom += math.exp((-1/(2*(float(Sigma**2))))*(float(X[0,i]-Mu[j]))**2)
        for j in xrange(0,k):
            Numer = math.exp((-1/(2*(float(Sigma**2))))*(float(X[0,i]-Mu[j]))**2)
            Expectations[i,j] = Numer / Denom
    if isdebug:
        print "***********"
        print u"隱藏變量E（Z）："
        print Expectations
# EM算法：步驟2，求最大化E[zij]的參數Mu
def m_step(k,N):
    global Expectations
    global X
    for j in xrange(0,k):
        Numer = 0
        Denom = 0
        for i in xrange(0,N):
            Numer += Expectations[i,j]*X[0,i]
            Denom +=Expectations[i,j]
        Mu[j] = Numer / Denom 
# 算法迭代iter_num次。或達到精度Epsilon停止迭代
def run(Sigma,Mu1,Mu2,k,N,iter_num,Epsilon):
    ini_data(Sigma,Mu1,Mu2,k,N)
    print u"初始<u1,u2>:", Mu
    for i in range(iter_num):
        Old_Mu = copy.deepcopy(Mu)
        e_step(Sigma,k,N)
        m_step(k,N)
        print i,Mu
        if sum(abs(Mu-Old_Mu)) < Epsilon:
            break
if __name__ == '__main__':
   run(6,40,20,2,1000,1000,0.0001)
   plt.hist(X[0,:],50)
   plt.show()

       本代碼用於模擬k=2個正態分布的均值預計。當中ini_data(Sigma,Mu1,Mu2,k,N)函數用於生成訓練樣本，此訓練樣本時從兩個高斯分布中隨機生成的，當中高斯分布a均值Mu1=40、均方差Sigma=6。高斯分布b均值Mu2=20、均方差Sigma=6，生成的樣本分布例如以下圖所看到的。

因為本問題中實現無法直接沖樣本數據中獲知兩個高斯分布參數，因此須要使用EM算法估算出詳細Mu1、Mu2取值。

圖 1  樣本數據分布

      在圖1的樣本數據下，在第11步時，迭代終止。EM預計結果為：

                                            Mu=[ 40.55261688  19.34252468]

附：

                                                    極大似然預計

參考文獻：機器學習TomM.Mitchell P.137