EM算法求高斯混合模型參数预计——Python实现

最新推荐文章于 2020-11-19 14:46:32 发布

weixin_33938733

最新推荐文章于 2020-11-19 14:46:32 发布

阅读量317

点赞数

文章标签： 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取值。