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算法估算出具体Mu一、Mu2取值。

图 1  样本数据分布

      在图1的样本数据下，在第11步时，迭代终止，EM估计结果为：

                                            Mu=[ 40.55261688  19.34252468]

附：

                                                    极大似然估计

参考文献：机器学习TomM.Mitchell P.137