[机器学习]如何应付EM算法考试

一、机器学习中的参数估计问题

    在前面的博文中,如“简单易学的机器学习算法——Logistic回归”中,采用了极大似然函数对其模型中的参数进行估计,简单来讲即对于一系列样本Logistic回归问题属于监督型学习问题,样本中含有训练的特征以及标签,在Logistic回归的参数求解中,通过构造样本属于类别和类别的概率:



这样便能得到Logistic回归的属于不同类别的概率函数:


此时,使用极大似然估计便能够估计出模型中的参数。但是,如果此时的标签是未知的,称为隐变量,如无监督的学习问题,典型的如K-Means聚类算法,此时不能直接通过极大似然估计估计出模型中的参数。

二、EM算法简介

    在上述存在隐变量的问题中,不能直接通过极大似然估计求出模型中的参数,EM算法是一种解决存在隐含变量优化问题的有效方法。EM算法是期望极大(Expectation Maximization)算法的简称,EM算法是一种迭代型的算法,在每一次的迭代过程中,主要分为两步:即求期望(Expectation)步骤和最大化(Maximization)步骤。

三、EM算法推导的准备

1、凸函数

    设是定义在实数域上的函数,如果对于任意的实数,都有

那么是凸函数。若不是单个实数,而是由实数组成的向量,此时,如果函数Hesse矩阵是半正定的,即


那么是凸函数。特别地,如果或者,那么称为严格凸函数。

2、Jensen不等式

    如果函数是凸函数,是随机变量,那么


特别地,如果函数是严格凸函数,那么当且仅当

即随机变量是常量。


(图片来自参考文章1)

注:若函数是凹函数,上述的符号相反。

3、数学期望

3.1随机变量的期望

   设离散型随机变量的概率分布为:


其中,,如果绝对收敛,则称的数学期望,记为,即:


   若连续型随机变量的概率密度函数为,则数学期望为:


3.2随机变量函数的数学期望

   设是随机变量的函数,即,若是离散型随机变量,概率分布为:


则:


   若是连续型随机变量,概率密度函数为,则


四、EM算法的求解过程

    假设表示观测变量,表示潜变量,则此时即为完全数据,的似然函数为,其中,为需要估计的参数,那么对于完全数据,的似然函数为
    构建好似然函数,对于给定的观测数据,为了估计参数,我们可以使用极大似然估计的方法对其进行估计。因为变量是未知的,我们只能对的似然函数为进行极大似然估计,即需要极大化:

上述式子中无法直接对求极大值,因为在函数中存在隐变量,即未知变量。若此时,我们能够确定隐变量的值,便能够求出的极大值,可以用过不断的修改隐变量的值,得到新的的极大值。这便是EM算法的思路。通过迭代的方式求出参数
    首先我们需要对参数赋初值,进行迭代运算,假设第次迭代后参数的值为,此时的log似然函数为,即:

在上式中,第二行到第三行使用到了Jensen不等式,由于log函数是凹函数,由Jensen不等式得到:



表示的是的期望,其中,表示的是隐变量满足的某种分布。这样,上式的值取决于的概率。在迭代的过程中,调整这两个概率,使得下界不断的上升,这样就能求得的极大值。注意,当等式成立时,说明此时已经等价于。由Jensen不等式可知,等式成立的条件是随机变量是常数,即:

已知:

所以:

则:

至此,我们得出了隐变量满足的分布的形式。这就是EM算法中的E步。在确定了后,调整参数使得取得极大,这便是M步。EM算法的步骤为:
  1. 初始化参数,开始迭代;
  2. E步:假设为第次迭代参数的估计值,则在第次迭代中,计算
  3. M步:求使极大化的,确定次的参数的估计值

五、EM算法的收敛性保证

迭代的过程能否保证最后找到的就是最大的似然函数值呢?即需要证明在整个迭代的过程中,极大似然估计是单调增加的。假定是EM算法的第次和第次迭代后的结果,选定,进行迭代:
  1. E步:
  2. M步:
固定,将看成变量:

上式中,第一个大于等于是因为:

六、利用EM算法参数求解实例

    假设有有一批数据分别是由两个正态分布:


产生,其中,未知,。但是不知道具体的是第产生,即可以使用表示。这是一个典型的涉及到隐藏变量的例子,隐藏变量为。可以使用EM算法对参数进行估计。

  1. 首先是初始化
  2. E步:,即求数据是由第个分布产生的概率:
  3. M步:,即计算最大的期望值。然而我们要求的参数是均值,可以通过如下的方式估计:

一个例子

高斯混合模型(Gaussian Mixed Model)指的是多个高斯分布函数的线性组合,理论上GMM可以拟合出任意类型的分布,通常用于解决同一集合下的数据包含多个不同的分布的情况(或者是同一类分布但参数不一样,或者是不同类型的分布,比如正态分布和伯努利分布)。

如图1,图中的点在我们看来明显分成两个聚类。这两个聚类中的点分别通过两个不同的正态分布随机生成而来。但是如果没有GMM,那么只能用一个的二维高斯分布来描述图1中的数据。图1中的椭圆即为二倍标准差的正态分布椭圆。这显然不太合理,毕竟肉眼一看就觉得应该把它们分成两类。

图1
图1

这时候就可以使用GMM了!如图2,数据在平面上的空间分布和图1一样,这时使用两个二维高斯分布来描述图2中的数据,分别记为N(μ1,Σ1)N(μ1,Σ1)N(μ2,Σ2)N(μ2,Σ2). 图中的两个椭圆分别是这两个高斯分布的二倍标准差椭圆。可以看到使用两个二维高斯分布来描述图中的数据显然更合理。实际上图中的两个聚类的中的点是通过两个不同的正态分布随机生成而来。如果将两个二维高斯分布N(μ1,Σ1)N(μ1,Σ1)N(μ2,Σ2)N(μ2,Σ2)合成一个二维的分布,那么就可以用合成后的分布来描述图2中的所有点。最直观的方法就是对这两个二维高斯分布做线性组合,用线性组合后的分布来描述整个集合中的数据。这就是高斯混合模型(GMM)。

图2
图2

高斯混合模型(GMM)

设有随机变量XX,则混合高斯模型可以用下式表示: 

p(x)=k=1KπkN(x|μk,Σk)p(x)=∑k=1KπkN(x|μk,Σk)

其中N(x|μk,Σk)N(x|μk,Σk)称为混合模型中的第kk分量(component)。如前面图2中的例子,有两个聚类,可以用两个二维高斯分布来表示,那么分量数K=2K=2πkπk混合系数(mixture coefficient),且满足: 

k=1Kπk=1∑k=1Kπk=1

0πk10≤πk≤1

可以看到πkπk相当于每个分量N(x|μk,Σk)N(x|μk,Σk)的权重。

GMM的应用

GMM常用于聚类。如果要从 GMM 的分布中随机地取一个点的话,实际上可以分为两步:首先随机地在这 K 个 Component 之中选一个,每个 Component 被选中的概率实际上就是它的系数πkπk ,选中 Component 之后,再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布,转化为已知的问题。

将GMM用于聚类时,假设数据服从混合高斯分布(Mixture Gaussian Distribution),那么只要根据数据推出 GMM 的概率分布来就可以了;然后 GMM 的 K 个 Component 实际上对应KK个 cluster 。根据数据来推算概率密度通常被称作 density estimation 。特别地,当我已知(或假定)概率密度函数的形式,而要估计其中的参数的过程被称作『参数估计』。

例如图2的例子,很明显有两个聚类,可以定义K=2K=2. 那么对应的GMM形式如下: 

p(x)=π1N(x|μ1,Σ1)+π2N(x|μ2,Σ2)p(x)=π1N(x|μ1,Σ1)+π2N(x|μ2,Σ2)

上式中未知的参数有六个:(π1,μ1,Σ1;π2,μ2,Σ2)(π1,μ1,Σ1;π2,μ2,Σ2). 之前提到GMM聚类时分为两步,第一步是随机地在这KK个分量中选一个,每个分量被选中的概率即为混合系数πkπk. 可以设定π1=π2=0.5π1=π2=0.5,表示每个分量被选中的概率是0.5,即从中抽出一个点,这个点属于第一类的概率和第二类的概率各占一半。但实际应用中事先指定πkπk的值是很笨的做法,当问题一般化后,会出现一个问题:当从图2中的集合随机选取一个点,怎么知道这个点是来自N(x|μ1,Σ1)N(x|μ1,Σ1)还是N(x|μ2,Σ2)N(x|μ2,Σ2)呢?换言之怎么根据数据自动确定π1π1π2π2的值?这就是GMM参数估计的问题。要解决这个问题,可以使用EM算法。通过EM算法,我们可以迭代计算出GMM中的参数:(πk,xk,Σk)(πk,xk,Σk).

GMM参数估计过程

GMM的贝叶斯理解

在介绍GMM参数估计之前,我们先改写GMM的形式,改写之后的GMM模型可以方便地使用EM估计参数。GMM的原始形式如下: 

p(x)=k=1KπkN(x|μk,Σk)(1)(1)p(x)=∑k=1KπkN(x|μk,Σk)

前面提到πkπk可以看成是第kk类被选中的概率。我们引入一个新的KK维随机变量zzzk(1kK)zk(1≤k≤K)只能取0或1两个值;zk=1zk=1表示第kk类被选中的概率,即:p(zk=1)=πkp(zk=1)=πk;如果zk=0zk=0表示第kk类没有被选中的概率。更数学化一点,zkzk要满足以下两个条件: 

zk{0,1}zk∈{0,1}

Kzk=1∑Kzk=1

例如图2中的例子,有两类,则zz的维数是2. 如果从第一类中取出一个点,则z=(1,0)z=(1,0);,如果从第二类中取出一个点,则z=(0,1)z=(0,1).

zk=1zk=1的概率就是πkπk,假设zkzk之间是独立同分布的(iid),我们可以写出zz的联合概率分布形式: 

p(z)=p(z1)p(z2)...p(zK)=k=1Kπzkk(2)(2)p(z)=p(z1)p(z2)...p(zK)=∏k=1Kπkzk

因为zkzk只能取0或1,且zz中只能有一个zkzk为1而其它zj(jk)zj(j≠k)全为0,所以上式是成立的。

图2中的数据可以分为两类,显然,每一類中的数据都是服从正态分布的。这个叙述可以用条件概率来表示: 

p(x|zk=1)=N(x|μk,Σk)p(x|zk=1)=N(x|μk,Σk)

即第kk类中的数据服从正态分布。进而上式有可以写成如下形式: 
p(x|z)=k=1KN(x|μk,Σk)zk(3)(3)p(x|z)=∏k=1KN(x|μk,Σk)zk

上面分别给出了p(z)p(z)p(x|z)p(x|z)的形式,根据条件概率公式,可以求出p(x)p(x)的形式: 

p(x)=zp(z)p(x|z)=i=1KπkN(x|μk,Σk) (zk=01)(4)(4)p(x)=∑zp(z)p(x|z)=∑i=1KπkN(x|μk,Σk) (zk=0的项为1,省略)

可以看到GMM模型的(1)式与(4)式有一样的形式,且(4)式中引入了一个新的变量zz,通常称为隐含变量(latent variable)。对于图2中的数据,『隐含』的意义是:我们知道数据可以分成两类,但是随机抽取一个数据点,我们不知道这个数据点属于第一类还是第二类,它的归属我们观察不到,因此引入一个隐含变量zz来描述这个现象。

注意到在贝叶斯的思想下,p(z)p(z)是先验概率, p(x|z)p(x|z)是似然概率,很自然我们会想到求出后验概率p(z|x)p(z|x): 

γ(zk)=p(zk=1|x)=p(zk=1)p(x|zk=1)p(x,zk=1)=p(zk=1)p(x|zk=1)Kj=1p(zj=1)p(x|zj=1) ()=πkN(x|μk,Σk)Kj=1πjN(x|μj,Σj) ((3)(4))(5)(5)γ(zk)=p(zk=1|x)=p(zk=1)p(x|zk=1)p(x,zk=1)=p(zk=1)p(x|zk=1)∑j=1Kp(zj=1)p(x|zj=1) (全概率公式)=πkN(x|μk,Σk)∑j=1KπjN(x|μj,Σj) (结合(3)(4))

上式中我们定义符号γ(zk)γ(zk)来表示来表示第kk个分量的后验概率。在贝叶斯的观点下,πkπk可视为zk=1zk=1的先验概率。

上述内容改写了GMM的形式,并引入了隐含变量zz和已知xx后的的后验概率γ(zk)γ(zk),这样做是为了方便使用EM算法来估计GMM的参数。

EM算法估计GMM参数

EM算法(Expectation-Maximization algorithm)分两步,第一步先求出要估计参数的粗略值,第二步使用第一步的值最大化似然函数。因此要先求出GMM的似然函数。

假设x={x1,x2,...,xN}x={x1,x2,...,xN},对于图2,xx是图中所有点(每个点有在二维平面上有两个坐标,是二维向量,因此x1,x2x1,x2等都用粗体表示)。GMM的概率模型如(1)式所示。GMM模型中有三个参数需要估计,分别是ππμμΣΣ. 将(1)式稍微改写一下: 

p(x|π,μ,Σ)=k=1KπkN(x|μk,Σk)(6)(6)p(x|π,μ,Σ)=∑k=1KπkN(x|μk,Σk)

为了估计这三个参数,需要分别求解出这三个参数的最大似然函数。先求解μkμk的最大似然函数。对(6)式取对数后再对μkμk求导并令导数为0即得到最大似然函数。 

0=n=1NπkN(xn|μk,Σk)jπjN(xn|μj,Σj)Σk(xnμk)(7)(7)0=−∑n=1NπkN(xn|μk,Σk)∑jπjN(xn|μj,Σj)Σk(xn−μk)

注意到上式中分数的一项的形式正好是(5)式后验概率的形式。两边同乘Σ1kΣk−1,重新整理可以得到: 

μk=1Nkn=1Nγ(znk)xn(8)(8)μk=1Nk∑n=1Nγ(znk)xn

其中: 

Nk=n=1Nγ(znk)(9)(9)Nk=∑n=1Nγ(znk)

(8)式和(9)式中,NN表示点的数量。γ(znk)γ(znk)表示点nnxnxn)属于聚类kk的后验概率。则NkNk可以表示属于第kk个聚类的点的数量。那么μkμk表示所有点的加权平均,每个点的权值是Nn=1γ(znk)∑n=1Nγ(znk),跟第kk个聚类有关。

同理求ΣkΣk的最大似然函数,可以得到: 

Σk=1Nkn=1Nγ(znk)(xnμk)(xnμk)T(10)(10)Σk=1Nk∑n=1Nγ(znk)(xn−μk)(xn−μk)T

最后剩下πkπk的最大似然函数。注意到πkπk有限制条件Kk=1πk=1∑k=1Kπk=1,因此我们需要加入拉格朗日算子: 

lnp(x|π,μ,Σ)+λ(k=1Kπk1)ln⁡p(x|π,μ,Σ)+λ(∑k=1Kπk−1)

求上式关于πkπk的最大似然函数,得到: 

0=n=1NN(xn|μk,Σk)jπjN(xn|μj,Σj)+λ(11)(11)0=∑n=1NN(xn|μk,Σk)∑jπjN(xn|μj,Σj)+λ

上式两边同乘πkπk,可以得到λ=Nλ=−N,进而可以得到πkπk更简洁的表达式: 

πk=NkN(12)(12)πk=NkN

EM算法估计GMM参数即最大化(8),(10)和(12)。需要用到(5),(8),(10)和(12)四个公式。我们先指定ππμμΣΣ的初始值,带入(5)中计算出γ(znk)γ(znk),然后再将γ(znk)γ(znk)带入(8),(10)和(12),求得πkπkμkμkΣkΣk;接着用求得的πkπkμkμkΣkΣk再带入(5)得到新的γ(znk)γ(znk),再将更新后的γ(znk)γ(znk)带入(8),(10)和(12),如此往复,直到算法收敛。

EM算法

  1. 定义分量数目KK,对每个分量kk设置πkπkμkμkΣkΣk的初始值,然后计算(6)式的对数似然函数。
  2. E step 
    根据当前的πkπkμkμkΣkΣk计算后验概率γ(znk)γ(znk) 
    γ(znk)=πkN(xn|μn,Σn)Kj=1πjN(xn|μj,Σj)γ(znk)=πkN(xn|μn,Σn)∑j=1KπjN(xn|μj,Σj)
  3. M step 
    根据E step中计算的γ(znk)γ(znk)再计算新的πkπkμkμkΣkΣk 
    μnewkΣnewkπnewk=1Nkn=1Nγ(znk)xn=1Nkn=1Nγ(znk)(xnμnewk)(xnμnewk)T=NkNμknew=1Nk∑n=1Nγ(znk)xnΣknew=1Nk∑n=1Nγ(znk)(xn−μknew)(xn−μknew)Tπknew=NkN

    其中: 
    Nk=n=1Nγ(znk)Nk=∑n=1Nγ(znk)
  4. 计算(6)式的对数似然函数 
    lnp(x|π,μ,Σ)=n=1Nln{k=1KπkN(xk|μk,Σk)}ln⁡p(x|π,μ,Σ)=∑n=1Nln⁡{∑k=1KπkN(xk|μk,Σk)}
  5. 检查参数是否收敛或对数似然函数是否收敛,若不收敛,则返回第2步。

Reference

  1. 漫谈 Clustering (3): Gaussian Mixture Model
  2. Draw Gaussian distribution ellipse
  3. Pang-Ning Tan 等, 数据挖掘导论(英文版), 机械工业出版社, 2010
  4. Christopher M. Bishop etc., Pattern Recognition and Machine Learning, Springer, 2006

Python代码

[python] view plain copy
  1. #coding:UTF-8  
  2. ''''' 
  3. Created on 2015年6月7日 
  4.  
  5. @author: zhaozhiyong 
  6. '''  
  7. from __future__ import division  
  8. from numpy import *  
  9. import math as mt  
  10. #首先生成一些用于测试的样本  
  11. #指定两个高斯分布的参数,这两个高斯分布的方差相同  
  12. sigma = 6  
  13. miu_1 = 40  
  14. miu_2 = 20  
  15.   
  16. #随机均匀选择两个高斯分布,用于生成样本值  
  17. N = 1000  
  18. X = zeros((1, N))  
  19. for i in xrange(N):  
  20.     if random.random() > 0.5:#使用的是numpy模块中的random  
  21.         X[0, i] = random.randn() * sigma + miu_1  
  22.     else:  
  23.         X[0, i] = random.randn() * sigma + miu_2  
  24.   
  25. #上述步骤已经生成样本  
  26. #对生成的样本,使用EM算法计算其均值miu  
  27.   
  28. #取miu的初始值  
  29. k = 2  
  30. miu = random.random((1, k))  
  31. #miu = mat([40.0, 20.0])  
  32. Expectations = zeros((N, k))  
  33.   
  34. for step in xrange(1000):#设置迭代次数  
  35.     #步骤1,计算期望  
  36.     for i in xrange(N):  
  37.         #计算分母  
  38.         denominator = 0  
  39.         for j in xrange(k):  
  40.             denominator = denominator + mt.exp(-1 / (2 * sigma ** 2) * (X[0, i] - miu[0, j]) ** 2)  
  41.           
  42.         #计算分子  
  43.         for j in xrange(k):  
  44.             numerator = mt.exp(-1 / (2 * sigma ** 2) * (X[0, i] - miu[0, j]) ** 2)  
  45.             Expectations[i, j] = numerator / denominator  
  46.       
  47.     #步骤2,求期望的最大  
  48.     #oldMiu = miu  
  49.     oldMiu = zeros((1, k))  
  50.     for j in xrange(k):  
  51.         oldMiu[0, j] = miu[0, j]  
  52.         numerator = 0  
  53.         denominator = 0  
  54.         for i in xrange(N):  
  55.             numerator = numerator + Expectations[i, j] * X[0, i]  
  56.             denominator = denominator + Expectations[i, j]  
  57.         miu[0, j] = numerator / denominator  
  58.           
  59.       
  60.     #判断是否满足要求  
  61.     epsilon = 0.0001  
  62.     if sum(abs(miu - oldMiu)) < epsilon:  
  63.         break  
  64.       
  65.     print step  
  66.     print miu  
  67.       
  68. print miu  

最终结果

[[ 40.49487592  19.96497512]]


参考文章:

1、(EM算法)The EM Algorithm (http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html)

2、数学期望(http://wenku.baidu.com/view/915a9c1ec5da50e2524d7f08.html?re=view)


阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页