Expectation Maximization(EM)算法

概率模型有时既含有观测变量（observable variable），又含有隐变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计或贝叶斯估计来计算模型参数。但是，当模型含有隐变量时，就不能简单地使用以上估计方法，而EM算法就是针对含有隐变量的概率模型参数的极大似然估计法。

从掷硬币实验说起

给定两枚硬币A和B，它们出现正面的概率分别为θA和θB(未知)，我们的目标是通过重复5次以下实验来估计θ=(θA,θB)。

选择一枚硬币，并且知道具体是哪一枚，然后用选中的硬币掷十次，记录结果。因此，总共进行了50次掷硬币的实验。如下图所示

Expectation Maximization(EM)算法

现在，引入两个随机变量X={X1,X2,X3,X4,X5}和Z={Z1,Z2,Z3,Z4,Z5}，其中Xi∈{0,1,⋯,10}表示第i次实验硬币正面朝上的次数，Zi∈{A,B}表示第i次实验所用的硬币。上述问题为完全数据的参数估计问题，可通过正面出现的比例进行估计

θ A^= # o f h e a d s u s i n g c o i n A t o t a l # o f f l i p s u s i n g c o i n A θ B^= # o f h e a d s u s i n g c o i n B t o t a l # o f f l i p s u s i n g c o i n B

事实上，上述的估计方式就是统计学上的最大似然估计。

现在考虑一个更有挑战性的参数估计问题。只给定硬币出现正面的次数，而不给定是由哪一枚硬币掷出的，即Z为隐含变量。因此，该问题就转化为不完全数据的参数估计问题。此时，由于不知道具体是哪一枚硬币，所以无法通过直接计算硬币出现正面的次数来估计θ。当概率模型存在隐变量时，不能简单地使用极大似然估计，需要采用下文的EM算法来计算模型参数。

凸函数

定义1 定义在区间I=[a,b]上的实函数f，如果对于∀x1,x2∈I,λ∈[0,1]，函数f满足以下不等式

f (λ x 1 + (1 - λ) x 2) \leq λ f (x 1) + (1 - λ) f (x 2) (1)

则称函数f为凸函数。其几何解释如下图所示
Expectation Maximization(EM)算法

定义2 如果函数−f是凸函数，则函数f是凹函数。

定理1 如果函数f在区间[a,b]上二阶可导，并且f′′(x)≥0，则f在区间[a,b]上为凸函数。

定理2(Jensen inequality) f为区间[a,b]上的凸函数，若x1,x2,⋯,xn∈I，并且λ1,λ2,⋯,λn≥0,∑ni=1=1，则以下不等式成立

f (\sum i = 1 n λ i x i) \leq \sum i = 1 n λ i f (x i) (2)

证明： 使用数学归纳法证明

当n=1时，(2)式显然成立，并取得等号。
当n=2时，即为凸函数的定义。
假设当n=n时(2)式成立，则

f (\sum i = 1 n + 1 λ i x i) = f (λ n + 1 x n + 1 + \sum i = 1 n λ i x i) = f (λ n + 1 x n + 1 + (1 - λ n + 1) 1 1 - λ n + 1 \sum i = 1 n λ i x i) \leq λ n + 1 f (x n + 1) + (1 - λ n + 1) f (1 1 - λ n + 1 \sum i = 1 n λ i x i)) = λ n + 1 f (x n + 1) + (1 - λ n + 1) f (\sum i = 1 n λ i 1 - λ n + 1 x i)) \leq λ n + 1 f (x n + 1) + (1 - λ n + 1) \sum i = 1 n λ i 1 - λ n + 1 f (x i) = λ n + 1 f (x n + 1) + \sum i = 1 n λ i f (x i) = \sum i = 1 n + 1 λ i f (x i)

得证。

因为λ1,λ2,⋯,λn≥0,∑ni=1=1，所以当X为随机变量时，Jensen inequality可以写成

E f (X) \geq f (E X) (3)

图形表示为

Expectation Maximization(EM)算法

EM算法

对于含有隐含变量的数据，EM算法是一个用于计算最大概率的高效的迭代过程，我们希望通过估计得到的模型参数可以使得观测到的数据出现的可能性最大。

EM算法的导出

对于一个含有隐含变量Z的概率模型，假设有m个样本{x(1),x(2),⋯,x(m)}，目标是极大化观测数据（不完全数据）X关于参数θ的对数似然函数，即极大化

l (θ) = \sum i = 1 m log p (x; θ) = \sum i = 1 m log \sum Z p (x, z; θ) (4)

极大化的主要困难在于上式中存在隐变量。最大化的l(θ)很难求得，但是我们可以构造l(θ)的下边界(E-step)，然后优化该下边界(M-step)，不断重复以上过程。

令Qi表示z的分布，满足条件∑zQi(z)=1,Qi(z)≥0

\sum i log p (x (i); θ) = \sum i log \sum z (i) p (x (i), z (i); θ) = \sum i log \sum z (i) Q i (z (i)) p (x (i), z (i); θ) Q i (z (i)) \geq \sum i \sum z (i) Q i (z (i)) log p (x (i), z (i); θ) Q i (z (i)) (5)

最后一步是由Jensen inequality得到的。并且，∑z(i)Qi(z(i))logp(x(i),z(i);θ)Qi(z(i))为p(x(i),z(i);θ)Qi(z(i))关于z(i)的期望。

对于任意的分布Qi，(5)式给出了l(θ)的下边界。现在的问题是如何选择Qi。

最直接的想法是，对于特定的θ，使得下边界尽可能地靠近l(θ)，即选择合适的Qi，使得(5)式的等号成立。由Jensen inequality可知，若要使得等号成立，需满足以下条件

p (x (i), z (i); θ) Q i (z (i)) = c

其中c不依赖于z(i)。因此，只要使得

Q i (z (i)) \propto p (x (i), z (i); θ)

因为Qi需满足∑zQi(z(i))=1，因此可以取

Q i (z (i)) = p (x (i), z (i); θ) \sum z p (x (i), z; θ) = p (x (i), z (i); θ) p (x (i); θ) = p (z (i) | x (i); θ) (6) (7)

(6)式等号右边对z(i)求和可以保证为1。

现在，通过将Qi设置为(7)式，我们得到了(5)式的最大化下边界，这就是E-step。在M-step中，求得一个新的θ，使得(5)式最大化。即

R e p e a t u n t i l c o n v e r g e n c e {(E - s t e p) F o r e a c h i, s e t Q i (z (i)) : = p (z (i) | x (i); θ) . (M - s t e p) S e t θ : = arg max θ \sum i \sum z (i) Q i (z (i)) log p (x (i), z (i); θ) Q i (z (i))}

下图为迭代过程的图形化表示

Expectation Maximization(EM)算法

收敛性

假设参数θ(t)和θ(t+1)是由两个连续的迭代过程求得的，下面来证明EM算法的每一次迭代都会改善似然函数，即l(θ(t))≤l(θ(t+1))。这个过程依赖Qi的选择，当给定参数θ(t)时，设置Qi为Q(t)i:=p(z(i)|x(i);θ(t))，由此可通过(5)式得到

l (θ (t)) = \sum i \sum z (i) Q (t) i (z (i)) log p (x (i), z (i); θ (t)) Q i (z (i)) (8)

同样由式(5)可得

l (θ (t + 1)) \geq \sum i \sum z (i) Q (t) i (z (i)) log p (x (i), z (i); θ (t + 1)) Q (t) i (z (i)) \geq \sum i \sum z (i) Q (t) i (z (i)) log p (x (i), z (i); θ (t)) Q (t) i (z (i)) = l (θ (t)) (9)

式(9)是因为θ(t+1)是由下式迭代得到的

arg max θ \sum i \sum z (i) Q i (z (i)) log p (x (i), z (i); θ) Q i (z (i))

综上所述，l(θ)单调递增。可以将在连续的两次迭代中，l(θ)值的变化小于某一阈值作为停止迭代的条件。

note： EM算法只能保证参数估计序列收敛到对数似然函数序列的极大值点，不能保证收敛到最大值点。所以在应用中，初值的选择变得非常重要，常用的方法是选取几个不同的初值进行迭代，然后对得到的各个估计值加以比较，从中选择最好的。

再看掷硬币问题

其中一种迭代方式入下：

初始化参数θ(t)^=(θ(t)A^,θ(t)B^)，同时初始化硬币A和B出现的概率（等概率）。
假设初始化的参数是正确的，通过最大似然估计计算模型参数θ(t+1)^.
重复以上两个过程直至收敛。

Expectation Maximization(EM)算法

从理论公式推导高斯混合模型

随机变量X是由k个高斯分布混合而成，取得各个高斯分布的概率为ϕ1,ϕ2,⋯,ϕk，第i个高斯分布的均值为μi，方差为Σi。若观测到随机变量X的一系列样本{x(1),x(2),⋯,x(m)}，试估计参数ϕ,μ,Σ。

E-step

w (i) j = Q i (z (i) = j) = p (z (i) = j | x (i); ϕ, μ, Σ)

w(i)j表示第i个样本x(i)属于第j个分模型的概率。

M-step

将多项分布和高斯分布的参数代入，并求解ϕ,μ,Σ使得式(5)取最大值，即

\sum i = 1 m \sum z (i) Q i (z (i)) log p (x (i), z (i); ϕ, μ, Σ) Q i (z (i)) = \sum i = 1 m \sum j = 1 k Q i (z (i) = j) log p (x (i) | z (i) = j; μ, Σ) p (z (i) = j; ϕ) Q i (z (i) = j) = \sum i = 1 m \sum j = 1 k w (i) j log 1 (2 π) n / 2 | Σ j | 1 / 2 exp (- 12 (x (i)) - μ j) T Σ - 1 j (x (i)) - μ j) \cdot ϕ j w (i) j (10)

对高斯分布的均值μl求偏导

\nabla μ l \sum i = 1 m \sum j = 1 k w (i) j log 1 (2 π) n / 2 | Σ j | 1 / 2 exp (- 12 (x (i) - μ j) T Σ - 1 j (x (i) - μ j) \cdot ϕ j w (i) j = - \nabla μ l \sum i = 1 m \sum j = 1 k w (i) j 12 (x (i) - μ j) T Σ - 1 j (x (i) - μ j) = 12 \sum i = 1 m w (i) j \nabla μ l 2 (μ T l Σ - 1 l x (i) - μ T l Σ - 1 l μ l) = \sum i = 1 m w (i) j (Σ - 1 l x (i) - Σ - 1 l μ l)

令其等于0可得

μ l : = \sum m i = 1 w (i) l x (i) \sum m i = 1 w (i) l (11)

同理，对高斯分布的方差求偏导可得

Σ j : = \sum m i = 1 w (i) j (x (i) - μ j) (x (i) - μ j) T \sum m i = 1 w (i) j (12)

求多项分布的参数ϕ，删除(10)式中的常数项，得

\sum i = 1 m \sum j = 1 k w (i) j log ϕ j s . t . \sum j = 1 k ϕ j = 1

由于log(⋅)已经约束了ϕj>0，所以约束条件中没有该约束。

建立拉格朗日方程

L (ϕ) = \sum i = 1 m \sum j = 1 k w (i) j log ϕ j + β (\sum j = 1 k ϕ j - 1) (13)

对ϕj求偏导，并令其等于0

\partial L (ϕ) \partial ϕ j = \sum i = 1 m w (i) j ϕ j + β - β = \sum i = 1 m \sum j = 1 k w (i) j = \sum i = 1 m 1 = m ϕ j = 1 m \sum i = 1 m w (i) j (14)

至此，由EM算法求解出了高斯混合模型的所有参数(11)(12)(14)。

参考文献

[1] 李航. 统计学习方法[M]. 清华大学出版社, 2012
[2] Borman S. The expectation maximization algorithm-a short tutorial[J]. Submitted for publication, 2004: 1-9.
[3] Do C B. What is the expectation maximization algorithm?[J]. Nature Biotechnology, 2008, 26(8):897-9.
[4] Andrew Ng. CS229 Lecture notes: The EM algorithm.