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

    在前面的博文中，如“简单易学的机器学习算法——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算法的收敛性保证

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

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

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

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

一个例子

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

如图1，图中的点在我们看来明显分成两个聚类。这两个聚类中的点分别通过两个不同的正态分布随机生成而来。但是如果没有GMM，那么只能用一个的二维高斯分布来描述图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

高斯混合模型（GMM）

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

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

可以看到π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形式如下： 

上式中未知的参数有六个：(π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的原始形式如下： 

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

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

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

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

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

即第kk类中的数据服从正态分布。进而上式有可以写成如下形式： 

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

可以看到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)γ(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)式稍微改写一下： 

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

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

其中： 

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

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

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

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

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

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=1Nk∑n=1Nγ(znk)xn=1Nk∑n=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代码

最终结果

[[ 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)