EM算法详细例子及推导

期望最大化算法

binary_seach

December 3, 2012

Contents 1简介

2一个例子[6]

3数学基础[7]

3.1极大似然估计[1]. . . . . . . . . . . . . . . . . . . . . . . . . .

3.2期望最大化算法收敛性. . . . . . . . . . . . . . . . . . . . . . 4应用举例

4.1参数估计. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2聚类. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

122467777

1简介

期望最大化(ExpectationMaximization) 算法最初是由Ceppellini[2]等人1950年在讨论基因频率的估计的时候提出的。后来又被Hartley[3]和Baum[4]等人发展的更加广泛。目前引用的较多的是1977年Dempster[5]等人的工作。它主要用于从不完整的数据中计算最大似然估计。后来经过其他学者的发展，这个算法也被用于聚类等应用。

2一个例子[6]

本章节希望通过一个例子来解释期望最大化算法的来由以及合理性。考虑一个投掷硬币的实验：现在我们有两枚硬币A 和B ，这两枚硬币和普通的硬币不一样，他们投掷出正面的概率和投掷出反面的概率不一定相同。我们将A 和B 投掷出正面的概率分别记为θA 和θB 。我们现在独立地做5次试验：随机的从这两枚硬币中抽取1枚，投掷10次，统计出现正面的次数。那么我们就得到了如表格1的实验结果。

2A

3A

4B

5A 9847

Table 1:硬币投掷实验的结果

在这个实验中，我们记录两组随机变量X =(X 1, X 2, X 3, X 4, X 5) , Z =(Z 1, Z 2, Z 3, Z 4, Z 5) ，其中X i ∈{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}代表试验i 中出现正面的次数，Z i ∈{A, B }代表这次试验投掷的是硬币A 还是硬币B 。我们的目标是通过这个实验来估计θ=(θA , θB ) 的数值。这个实验中的参数估计就是有完整数据的参数估计，这个是因为我们不仅仅知道每次试验中投掷出正面的次数，我们还知道每次试验中投掷的是硬币A 还是B 。一个很简单也很直接的估计θ的方法如公式1所示。

用硬币A 投掷出正面的次数θˆA =用硬币A 投掷的次数

用硬币B 投掷出正面的次数θˆB =用硬币B 投掷的次数(1)

实际上这样的估计就是统计上的极大似然估计(maximumlikelihood estimation) 的结果。用P (X, Z |θ) 来表示X,Z 的联合概率分布（其中带有参

2

数θ），那么对于上面的实验，我们可以计算出他们出现我们观察到的结果即x 0=(5, 9, 8, 4, 7) , z 0=(B, A, A, B, A ) 的概率2。函数P (X =x (0), Z =z (0)|θ) 就叫做θ的似然函数。我们将它对θ求偏导并令偏导数为0，就可以得到如1的结果。

3P (Z =A ) 3(1−P (Z =A )) 2P (X =x (0), Z =z (0)|θ) =C 5

99×C 10θA (1−θA )

88θA (1−θA ) 2×C 10

77×C 10θA (1−θA ) 3

55×C 10θB (1−θB ) 5

44×C 10θB (1−θB ) 6(2)

我们将这个问题稍微改变一下，我们将我们所观察到的结果修改一下：我们现在只知道每次试验有几次投掷出正面，但是不知道每次试验投掷的是哪个硬币，也就是说我们只知道表1中第一列和第三列。这个时候我们就称Z 为隐藏变量(HiddenVariable) ，X 称为观察变量(ObservedVariable) 。这个时候再来估计参数θA 和θB ，就没有那么多数据可供使用了，这个时候的估计叫做不完整数据的参数估计。

如果我们这个时候有某种方法（比如，正确的猜到每次投掷硬币是A 还是B ），这样的话我们就可以将这个不完整的数据估计变为完整数据估计。

当然我们如果没有方法来获得更多的数据的话，那么下面提供了一种在这种不完整数据的情况下来估计参数θ的方法。我们用迭代的方式来进行：

(1)我们先赋给θ一个初始值，这个值不管是经验也好猜的也好，反正我

们给它一个初始值。在实际使用中往往这个初始值是有其他算法的结果给出的，当然随机给他分配一个符合定义域的值也可以。这里我们就给定θA =0. 7, θB =0. 4

(2)然后我们根据这个来判断或者猜测每次投掷更像是哪枚硬币投掷的结

果。比如对于试验1，如果投掷的是A ，那么出现5个正面的概率为C 105×0. 75×(1−0. 7) 5≈0. 1029；如果投掷的是B ，出现5个正面的概率为C 105×0. 45×(1−0. 4) 5≈0. 2007；基于试验1的试验结果，可以判断这个试验投掷的是硬币A 的概率为0.1029/(0.1029+0.2007)=0.3389，是B 的概率为0.2007/(0.1029+0.2007)=0.6611。因此这个结果更可能是投掷B 出现的结果。

(3)假设上一步猜测的结果为B,A,A,B,A ，那么根据这个猜测，可以像完

整数据的参数估计一样(公式2) 重新计算θ的值。

这样一次一次的迭代2-3步骤直到收敛，我们就得到了θ的估计。现在你可能有疑问，这个方法靠谱么？事实证明，它确实是靠谱的。

3

期望最大化算法就是在这个想法上改进的。它在估计每次投掷的硬币的时候，并不要确定住这次就是硬币A 或者B ，它计算出来这次投掷的硬币是A 的概率和是B 的概率；然后在用这个概率（或者叫做Z 的分布）来计算似然函数。期望最大化算法步骤总结如下：

E 步骤先利用旧的参数值θ′计算隐藏变量Z 的（条件）分布P θ′(Z =z j |X i =

x i ) ，然后计算logP θ(Z, X =x ) 1的期望

E (logP θ(Z, X =x )) =n ∑∑

i =1j P θ′(Z =z j |X i =x i ) logP θ(Z =z j , X i =x i )

(3)

其中θ是当前的θ值，而θ′是上一次迭代得到的θ值。公式3中已经只剩下θ一个变量了，θ′是一个确定的值，这个公式或者函数常常叫做Q 函数，用Q (θ,θ′) 来表示。

M 步骤极大化Q ，往往这一步是求导，得到由旧的θ值θ′来计算新的θ值的

公式∂Q=0(4)∂θ

总结一下，期望最大化算法就是先根据参数初值估计隐藏变量的分布，然后根据隐藏变量的分布来计算观察变量的似然函数，估计参数的值。前者通常称为E 步骤，后者称为M 步骤。

3数学基础[7]

首先来明确一下我们的目标：我们的目标是在观察变量X 和给定观察样本x 1, x 2,..., x n 的情况下，极大化对数似然函数

l (θ) =n ∑

i =1ln P (X i =x i ) (5)

其中只包含观察变量的概率密度函数

∑P (X i =x i ) =P (X i =x i , Z =z j )

j

1(6)这里因为参数θ的写法与条件概率的写法相同，因此将参数θ写到下标以更明确的表述

4

其中Z 为隐藏随机变量,{z j }是Z 的所有可能的取值。那么

l (θ) =

=n ∑i =1n ∑

i =1ln ∑j P (X i =x i , Z =z j ) P (X i =x i , Z =z j ) αj (7)ln ∑j αj

这里我们引入了一组参数（不要怕多，我们后面会处理掉它的）α，它满足∑∀可能的j, αj ∈(0, 1]和j αj =1到这里，先介绍一个凸函数的性质，或者∑叫做凸函数的定义。f (x ) 为凸函数，∀i =1, 2,..., n, λi ∈[0, 1], n

i =1λi =1，

对f (x ) 定义域中的任意n 个x 1, x 2,..., x n 有

n ∑∑f (λi x i ) ≤nλi f (x i )) 2

i =1i =1(8)

对于严格凸函数，上面的等号只有在x 1=x 2=... =x n 的时候成立。关于凸函数的其他性质不再赘述。对数函数是一个严格凸函数。因而我们可以有下面这个结果：

n ∑∑P (X i =x i , Z =z j ) l (θ) =ln αj αj i =1j (9)n ∑∑P (X i =x i , Z =z j ) ≥αj ln αj i =1j

现在我们根据等号成立的条件来确定αj 即

P (X i =x i , Z =z j ) =c (10)αj ∑其中c 是一个与j 无关的常数。因为j αj =1，稍作变换就可以得到

αj =P (X i =x i , Z =z j )

P (X i =x i ) (11)

现在来解释一下我们得到了什么。αj 就是Z =z j 在X =x i 下的条件概率或者后验概率。求α就是求隐藏随机变量Z 的条件分布。总结一下目前得到的公式就是

n ∑∑P (X i =x i , Z =z j ) l (θ) =αj ln (12)αj i =1j

直接就极大值比较难求，EM 算法就是按照下面这个过程来的。2它就是大名鼎鼎的琴生（Jensen ）不等式

5

(1)根据上一步的θ来计算α，即隐藏变量的条件分布

(2)极大化似然函数来得到当前的θ的估计

3.1极大似然估计[1]好吧，我觉得还是再说说极大似然估计吧。给定一个概率分布D ，假设其概率密度函数为f ，其中f 带有一组参数θ。为了估计这组参数θ，我们可以从这个分布中抽出一个具有n 个采样值的X 1, X 2,..., X n ，那么这个就是n 个（假设独立）同分布随机变量，他们分别有取值x 1, x 2,..., x n ，那么我们就可以计算出出现这样一组观察值的概率密度为

l (θ) =n ∏

i =1f (x i ) (13)

对于f 是离散的情况，就计算出现这一组观察值的概率。

l (θ) =n ∏

i =1P (X i =x i ) (14)

注意，这个函数中是含有参数θ的。θ的极大似然估计就是求让上面似然函数取极大值的时候的参数θ值。

一般来说，会将上面那个似然函数取自然对数，这样往往可以简化计算。记住，这样仅仅是为了简化计算3。取了自然对数之后的函数叫做对数似然函数。n ∑ln l (θ) =ln f (x i ) (15)

i =1

因为对数是一个严格单调递增的凹函数4，所以对似然函数取极大值与对对数似然函数取极大值是等价的。

取了对数之后还可以跟信息熵等概念联系起来

关于凸函数有很多种说法，上凸函数和下凸函数，凸函数和凹函数等等，这里指的是二阶导数大于（等于）0的一类函数，而凹函数是其相反数为凸函数的一类函数43

6

3.2期望最大化算法收敛性

n ∑∑

i =1j

(t ) αj 如何保证算法收敛呢？我们只用证明l (θ(t +1)) ≥l (θ(t ) ) 就可以了。l (θ(t +1)) =(t +1)αj ln ln P (X =x i , Z =z j ) (t +1)αj (t +1)≥n ∑∑

i =1j P (X =x i , Z =z j ) (t +1)αj

ln (t ) (16)≥n ∑∑

i =1j (t ) αj P (X =x i , Z =z j ) (t ) αj (t )

=l (θ(t ) )

其中第一个大于等于号是因为只有当α取值合适（琴生不等式等号成立条件）的时候才有等号成立，第二个大于等于号正是M 步骤的操作所致。这样我们就知道l (θ) 是随着迭代次数的增加越来越大的，收敛条件是值不再变化或者变化幅度很小。

4

4.1应用举例参数估计

很直接的应用就是参数估计，上面举的例子就是参数估计。

4.2聚类

但是如果估计的参数可以表明类别的话，比如某个参数表示某个样本是否属于某个集合。这样的话其实聚类问题也就可以归结为参数估计问题。References

[1]最大似然估计[Online].Available:http://zh.wikipedia.org/wiki/%E6%9C%80%E5%A4%A7%E4%BC%BC%E7%84%B6%E4%BC%B0%E8%AE%A1

[2]Ceppellini, r., Siniscalco, M. &Smith, C.A. Ann. Hum. Genet. 20, 97–115(1955).

[3]Hartley, H. Biometrics 14, 174–194(1958).

[4]Baum, L.E., Petrie, T., Soules, G. &Weiss, N. Ann. Math. Stat. 41, 164–171(1970).

7

[5]Dempster, A.P.; Laird, N.M.; Rubin, D.B. (1977)."Maximum Likelihood

from Incomplete Data via the EM Algorithm". Journal of the Royal Statis-tical Society. Series B (Methodological)39(1):1–38.JSTOR 2984875. MR 0501537.

[6]What is the expectation maximization algorithm? [Online].Avaiable:http:

//ai.stanford.edu/~chuongdo/papers/em_tutorial.pdf

[7]The EM Algorithm [Online].Available:http://www.cnblogs.com/

jerrylead/archive/2011/04/06/2006936.html

8