Opencv2.4.9源码分析——Expectation Maximization

一、原理

 

期望极大算法（EM，Expectation Maximization）是一种能够得到极大似然参数估计的迭代方法。虽然EM算法在不同的领域被更早的提出，但它是由Dempster等人于1977年被正式命名并给出解释。

对于一个除了拥有未知参数和可观测变量外，还包含隐含变量的统计模型来说，是无法用极大似然方法直接得到参数的，这是因为我们不能像极大似然法那样，同时对未知参数和隐含变量求导来求解似然函数。但该模型可以应用EM算法得到，它是通过可观测变量分别求解两个方程，即把第一个方程的解代入第二个方程中，再把第二个方程的解代入第一个方程中，依次类推直到收敛为止。虽然EM算法不能保证得到全局极值点，但可以通过一些方法得到改善，如选取不同的初始值，比较最终的结果。

下面我们就来详细推导EM算法。

给定一组相互独立且同分布的训练样本X＝{ x₁,x₂,…,x_n}，每个样本都是一个d维的向量，即每个样本都具有d个特征属性，x_i∈R^d，设p(X|θ)为被参数θ控制的概率密度函数，如X服从高斯正态分布，则θ表示均值μ和标准差σ，即X~N(μ, σ²)：

（1）

似然函数是一种关于统计模型中参数的函数，表示模型参数中的似然性，则参数θ在X下的似然函数为：

（2）

实际上为了便于计算，我们往往应用的是对数似然函数，即把式2写为：

（3）

该似然函数L(θ|X)反映了在基于可观测到的样本值X下，与真实统计模型相似程度的大小，因此我们有理由认为L(θ|X)的最大值所对应的参数θ^*就是真实模型的参数，即

（4）

这种参数估计的方法就称为极大似然法，而参数θ^*就是对L(θ|X)求一阶偏导为0时所对应的值。极大似然法对单一统计模型的参数估计简单实用，但它不适用对混合统计模型参数的估计。

设Z为某一样本数据，它除了有可观测的变量X＝{ x₁,x₂,…,x_n}外，还包含一个隐含变量Y＝{ y₁,y₂,…,y_n}，Y可以看成是X的输出或类别标签，y_i∈{1,…,K}，K表示样本类别的数量，这意味着Z的统计模型是混合统计模型。我们把Z=(X, Y)称为完整数据，而把X称为不完整数据。无论样本是否包含隐含变量，参数的估计都是基于所观测到的数据的，所以我们要对Z的统计模型参数θ进行估计，应用的是基于所观测到的不完整数据X的似然函数：

（5）

式中，p(x_i| θ)＝∑_jp(x_i,y_j|θ)表示观测到的变量x_i的概率分布为所有类标签y_i(i＝1,…,K)下x_i的概率之和。式5中包含了“和的对数（ln∑）”的形式，因此我们还是应用导数为0的方法求最大值就不那么容易了。

很明显，如果Y已知，则通过完整数据Z的似然函数——L(θ; X, Y)＝ln p(X, Y| θ)可以很容易的得到参数θ。L(θ; X, Y)其实是一个X和θ是常数，Y是随机变量的函数。

基于引入的隐含变量，EM算法给出了一种非常有效的处理极大似然参数估计的迭代方法。每次迭代都包括两个步骤：期望步骤（E-Step）和极大步骤（M-Step）。

在E-Step中，需要计算在给定X和当前参数估计下，关于Y的完整数据的对数似然函数ln p(X, Y| θ)的期望值，其定义为：

（6）

式中，θ^(t-1)为当前的参数估计，我们用该参数来得到期望值，而θ是新的参数，我们用它来优化提高Q值。关于式6，我们需要明确的是X和θ^(t-1)是常数，θ是我们希望调整的一般变量，而Y是基于分布函数p(Y| X, θ^(t-1))的离散随机变量，则由数学期望的定义，lnp(X, Y| θ)的期望值表示为：

（7）

p(Y| X, θ^(t-1))是依赖于可观测变量X和当前参数θ^(t-1)的、未观测变量Y的边缘分布函数，ϒ表示y_i的取值范围，即前面定义的{1,…,K}。

在M-Step中，最大化E-Step的期望值：

（8）

θ^(t)为本次迭代得到估计参数，它被用于下次迭代的E-Step中。那么第一次迭代中E-Step所使用的θ⁽⁰⁾为任意选取的。

EM算法的精髓就是用完整数据的对数似然函数的期望值不断迭代，最终得到不完整数据的最大似然估计。我们分析一下为什么会有这样的结果。

由联合分布概率可知：

（9）

则X的对数似然函数为：

（10）

对上式两边同乘以p(y_i | X, θ^(t-1))，则

（11）

然后对上式两边在y_i的取值范围ϒ内进行求和，因为

（12）

所以式11改写为：

（13）

上式第一项就是式7所定义的Q(θ; θ^(t-1))，而第二项我们定义为H(θ;θ^(t-1))，则

（14）

为什么我们仅仅通过迭代的方式不断找到能使Q(θ; θ^(t-1))达到最大值的θ^(t)就能保证此时的似然函数最大呢？即

（15）

式15表现的是EM算法的单调性。如果满足下列条件，式15就必然成立：

（16）

（17）

式16已经能够由EM算法中的M-Step得到保证（式8），所以关键就要证明式17成立。这里面我们要用到一个很重要的概念——Jensen不等式。

设f为定义在I＝[a, b]的一个实数域函数，如果满足下列条件，则f称为凸函数：

（18）

很显然，如果－f是凸函数，则f就是凹函数，ln(x)函数是凹函数。

如果f是定义在区间I上的凸函数，则满足：

（19）

式19的证明会用到式18的结论，下面我们就给出式19的证明：

（20）

我们称式19所表示的不等式为Jensen不等式，而当λ_i＝f(x_i)时，

（21）

因为ln(x)函数是凹函数，所以有

（22）

有了Jensen不等式的概念，我们再回到式17。要想使式17成立，只要证明下式即可：

（23）

而对于任意θ来说，下式成立，则式23就成立：

（24）

下面我们就来证明式24：

（25）

由于p(y_i | X, θ^(t-1))是概率分布，所以它满足Jensen不等式中对λ_i的条件（对比式12和式19中λ_i的条件），因此应用Jensen不等式，式25改写为：

（26）

式24成立，则式17就成立，最终式15就成立，即不断迭代得到的完整数据的似然函数的期望值可以保证收敛于局部极值。

EM算法的原理部分就介绍到这里，下面我们以高斯混合模型（GMM）为例，具体讲解如何由EM算法得到该模型的参数。

设高斯混合模型由K个高斯分布模型（即K个高斯成分）构成：

（27）

可以看出，高斯混合模型的参数θ不仅包含高斯参数：μ和∑，还包括每个高斯成分的权值因子α_k，因此参数θ为：

（28）

样本X＝{ x₁,x₂,…,x_n}，x_i∈R^d，隐含变量Y＝{ y₁,y₂,…,y_n}，y_i表示x_i属于K个高斯成分的哪一个。则式27可以进一步写为：

（29）

为了计算式6中的Q函数，我们还需要给出p(X, Y| θ)和p(Y| X, θ)。

（30）

式中，p(x_i| y_i,θ)表示y_i所指向的高斯成分的概率分布模型，p(y_i|θ)表示y_i所指向的高斯成分的权值因子α_k，因此p(X, Y|θ)又可写为

（31）

p(Y| X, θ)定义为

（32）

则式6中的Q函数表示为：

（33）

我们还需要对式33做进一步化简，此时我们用到了下面这个定理：

设p(Y)为联合概率分布：p(Y)=p(y₁, y₂,…,y_n)= p₁(y₁)p₂(y₂)…p_n(y_n)，F(Y)为一个线性函数，且F(Y)=f₁(y₁)+f₂(y₂)+…+f_n(y_n)，则

（34）

我们简单证明一下式34：

（35）

对式35中的累加和，我们仅拿出第一项来进行分析，即

（36）

很显然，式36中方括号内的部分就是p₁(y₁)的边缘概率分布，所以

（37）

我们把式37的结论应用到式35中累加和的其他项，则就得到了式34

（38）

式33和式34相比较，可以看出，式33中的log p(x_i, y_i|θ)就相当于式34中的f_i(y_i)，p(y_i| x_i,θ^(t-1))就相当于p_i(y_i)，所以Q函数又写为：

（39）

式39的第二个等式仅仅是交换了两个累加和的次序。我们用高斯模型分布代替p(x_i, y_i|θ)，并用混合高斯模型代替隐含变量y_i，则式39又写为：

（40）

我们通过一步步的推导终于得到了高斯混合模型EM算法中E-Step所需的Q函数，下面就进入M-Step，计算Q函数的参数θ的最大值。由式28可知，参数θ包括α_k，μ_k和∑_k。我们再对式40进行分解：