参考资料:
- 《统计学习方法》,李航著
- 知乎文章:【机器学习】EM——期望最大(非常详细)、人人都懂EM算法
文章目录
1 简述
EM算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法,含有隐变量的概率模型的数据表示为 P ( Y , Z ∣ θ ) P(Y,Z|\theta) P(Y,Z∣θ)。
这里,Y是观测变量的数据,Z是隐变量的数据, θ \theta θ是模型参数。
EM算法通过迭代求解观测数据的对数似然函数 L ( θ ) = l o g P ( Y ∣ θ ) L(\theta)=log P(Y|\theta) L(θ)=logP(Y∣θ)的极大化,实现极大似然估计
每次迭代包括两步:E步,求期望,即求
L
(
θ
)
=
l
o
g
P
(
Y
∣
θ
)
L(\theta)=log P(Y|\theta)
L(θ)=logP(Y∣θ)关于
P
(
Z
∣
Y
,
θ
(
i
)
)
P(Z|Y,\theta^{(i)})
P(Z∣Y,θ(i))的期望:
Q
(
θ
,
θ
(
i
)
)
=
∑
Z
l
o
g
P
(
Y
,
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
Q(\theta,\theta^{(i)})=\sum_{Z}log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})
Q(θ,θ(i))=∑ZlogP(Y,Z∣θ)P(Z∣Y,θ(i))
称为Q函数,这里 θ ( i ) \theta^{(i)} θ(i)是参数的现估计值;M步,求极大,即极大化Q函数得到参数的新估计值:
θ ( i + 1 ) = a r g m a x Q ( θ , θ ( i ) ) \theta^{(i+1)}=argmax Q(\theta,\theta^{(i)}) θ(i+1)=argmaxQ(θ,θ(i))
在构建具体的EM算法时,重要的是定义Q函数。
每次迭代中,EM算法通过极大化Q函数来增大对数似然函数 L ( θ ) L(\theta) L(θ)
EM算法在每次迭代后均提高观测数据的似然函数值,即 P ( Y ∣ θ ( i + 1 ) ) ≥ P ( Y ∣ θ ( i ) ) P(Y|\theta^{(i+1)})\ge P(Y|\theta^{(i)}) P(Y∣θ(i+1))≥P(Y∣θ(i))
在一般条件下EM算法是收敛的,但不能保证收敛到全局最优。
2 详述
EM 算法,全称 Expectation Maximization Algorithm。
EM 算法的核心思想非常简单,分为两步:Expection-Step 和 Maximization-Step。E-Step 主要通过观察数据和现有模型来估计参数,然后用这个估计的参数值来计算似然函数的期望值;而 M-Step 是寻找似然函数最大化时对应的参数。由于算法会保证在每次迭代之后似然函数都会增加,所以函数最终会收敛。
2.1 简单的例子
现在一个班里有 50 个男生和 50 个女生,且男女生分开。我们假定男生的身高服从正态分布: N ( μ 1 , σ 1 2 ) N(\mu_1,\sigma_1^2) N(μ1,σ12) ,女生的身高则服从另一个正态分布: N ( μ 2 , σ 2 2 ) N(\mu_2,\sigma_2^2) N(μ2,σ22)。这时候我们可以用极大似然法(MLE),分别通过这 50 个男生和 50 个女生的样本来估计这两个正态分布的参数。
但现在我们让情况复杂一点,就是这 50 个男生和 50 个女生混在一起了。我们拥有 100 个人的身高数据,却不知道这 100 个人每一个是男生还是女生。
这时候情况就有点尴尬,因为通常来说,我们只有知道了精确的男女身高的正态分布参数我们才能知道每一个人更有可能是男生还是女生。但从另一方面去考量,我们只有知道了每个人是男生还是女生才能尽可能准确地估计男女各自身高的正态分布的参数。
这个时候有人就想到我们必须从某一点开始,并用迭代的办法去解决这个问题:我们先设定男生身高和女生身高分布的几个参数(初始值),然后根据这些参数去判断每一个样本(人)是男生还是女生,之后根据标注后的样本再反过来重新估计参数。之后再多次重复这个过程,直至稳定。这个算法也就是 EM 算法。
2.2 基础知识
2.2.1 极大似然估计
前期准备
对于这个函数: p ( x ∣ θ ) p(x|\theta) p(x∣θ) 输入有两个:x表示某一个具体的数据; θ \theta θ表示模型的参数
如果 θ \theta θ是已知确定的, x是变量,这个函数叫做概率函数(probability function),它描述对于不同的样本点x,其出现概率是多少。
如果x是已知确定的, θ \theta θ是变量,这个函数叫做似然函数(likelihood function),它描述对于不同的模型参数,出现x这个样本点的概率是多少。
极大似然和EM(Expectation Maximization)算法,与其说是一种算法,不如说是一种解决问题的思想,解决一类问题的框架,和线性回归,逻辑回归,决策树等一些具体的算法不同,极大似然和EM算法更加抽象,是很多具体算法的基础。
极大似然估计的前提一定是要假设数据总体的分布,如果不知道数据分布,是无法使用极大似然估计的。
假设我们要调查学校学生的身高分布,假设学校学生的身高服从正态分布 N ( μ , σ 2 ) N(\mu,\sigma^2) N(μ,σ2)。
为方便统计,我们随机抽样,抽到了200个人,根据这200个人的身高估计均值 μ \mu μ和方差 σ 2 \sigma^2 σ2。
抽到某个学生的概率密度 p ( x ∣ θ ) p(x|\theta) p(x∣θ)服从高斯分布 N ( μ , σ 2 ) N(\mu,\sigma^2) N(μ,σ2),其中的未知参数是 θ = [ μ , θ ] T \theta=[\mu,\theta]^T θ=[μ,θ]T
同时抽到这200个学生的概率就是他们各自概率的乘积,用联合概率表示就是:
L ( θ ) = L ( x 1 , x 2 , . . . , x n ; θ ) = ∏ i = 1 n p ( x i , θ ) L(\theta)=L(x_1,x_2,...,x_n;\theta)=\prod_{i=1}^{n}p(x_i,\theta) L(θ)=L(x1,x2,...,xn;θ)=i=1∏np(xi,θ)
n为抽取的样本数,这个概率反映了,在概率密度函数的参数是 θ \theta θ 时,得到 X 这组样本的概率。上式中等式右侧只有 θ \theta θ是未知数,所以 L 是 θ \theta θ的函数。
换言之,这个函数反映的是在不同的参数 θ \theta θ取值下,取得当前这个样本集的可能性,也即是样本集的似然函数(likelihood function),记为 L ( θ ) L(\theta) L(θ)。
问题来了,怎么极大似然函数呢?
求 L ( θ ) L(\theta) L(θ) 对所有参数的偏导数,然后让这些偏导数为 0,假设有 n n n 个参数,就有 n n n 个方程组成的方程组,那么方程组的解就是似然函数的极值点了,从而得到对应的 L ( θ ) L(\theta) L(θ) 了。
总而言之,极大似然估计就是利用已知的样本结果信息,反推最具有可能(最大概率)导致这些样本结果出现的模型参数值。
2.2.2 Jensen不等式
如果函数f是凹函数,X是随机变量,则
E
[
f
(
X
)
]
≤
f
(
E
(
X
)
)
E[f(X)] \le f(E(X))
E[f(X)]≤f(E(X)),当f严格是凹函数时,则
E
[
f
(
X
)
]
<
f
(
E
(
X
)
)
E[f(X)] < f(E(X))
E[f(X)]<f(E(X)),凸函数反之。当
X
=
E
(
X
)
X=E(X)
X=E(X)时,等号成立。
l
2.2.3 期望
对于离散型的随机变量X的概率分布为
p
i
=
p
X
=
x
i
p_i=p{X=x_i}
pi=pX=xi,数学期望E(X)为:
E
(
x
)
=
∑
i
x
i
p
i
E(x)=\sum_{i}x_ip_i
E(x)=i∑xipi
p
i
p_i
pi是权值,满足两个条件
1
≥
p
i
≥
0
1\ge pi \ge 0
1≥pi≥0,
∑
i
p
i
=
1
\sum_i pi=1
∑ipi=1
若Y=g(x),若X是离散型随机变量,则:
E
(
Y
)
=
∑
i
g
(
x
i
)
p
i
E(Y)=\sum_ig(x_i)p_i
E(Y)=i∑g(xi)pi
2.3 EM算法
2.3.1 问题引入
如果我们现在已知的信息只有200个人的身高,至于具体每个人的性别未知。
这个时候,对于每一个样本或者你抽取到的人,就有两个问题需要估计了,一是这个人是男的还是女的,二是男生和女生对应的身高的正态分布的参数是多少。这两个问题是相互依赖的:
- 当我们知道了每个人是男生还是女生,我们可以很容易利用极大似然对男女各自的身高的分布进行估计。
- 反过来,当我们知道了男女身高的分布参数我们才能知道每一个人更有可能是男生还是女生。
但是现在我们既不知道每个学生是男生还是女生,也不知道男生和女生的身高分布。这就成了一个先有鸡还是先有蛋的问题了。鸡说,没有我,谁把你生出来的啊。蛋不服,说,没有我,你从哪蹦出来啊。为了解决这个你依赖我,我依赖你的循环依赖问题,总得有一方要先打破僵局,不管了,我先随便整一个值出来,看你怎么变,然后我再根据你的变化调整我的变化,然后如此迭代着不断互相推导,最终就会收敛到一个解(草原上的狼和羊,相生相克)。这就是EM算法的基本思想了。
EM的意思是“Expectation Maximization”,具体方法为:
- 先设定男生和女生的身高分布参数(初始值);
- 然后计算出每个人更可能属于第一个还是第二个正态分布中的(例如,这个人的身高是180,那很明显,他极大可能属于男生),这个是属于Expectation 一步;
- 我们已经大概地按上面的方法将这 200 个人分为男生和女生两部分,我们就可以根据之前说的极大似然估计分别对男生和女生的身高分布参数进行估计(这不变成了极大似然估计了吗?极大即为Maximization)这步称为 Maximization;
- 然后,当我们更新这两个分布的时候,每一个学生属于女生还是男生的概率又变了,那么我们就再需要调整E步;
- ……如此往复,直到参数基本不再发生变化或满足结束条件为止。
注意:
上述的学生属于男生还是女生我们称之为隐含参数,女生和男生的身高分布参数称为模型参数。
EM算法解决这个思路是使用启发式的迭代方法,既然我们无法直接求出模型分布参数,那么我们可以先猜想隐含参数(EM 算法的 E 步),接着基于观察数据和猜测的隐含参数一起来极大化对数似然,求解我们的模型参数(EM算法的M步)。由于我们之前的隐含参数是猜测的,所以此时得到的模型参数一般还不是我们想要的结果。我们基于当前得到的模型参数,继续猜测隐含参数(EM算法的 E 步),然后继续极大化对数似然,求解我们的模型参数(EM算法的M步)。以此类推,不断的迭代下去,直到模型分布参数基本无变化,算法收敛,找到合适的模型参数。
2.3.2 算法流程
后面的部分完全参考的《统计学习方法》。
记上图中的(9.10)式为式(1),后面会用到。
2.3.3 算法导出
上面叙述了EM算法。为什么EM算法能近似实现对观测数据的极大似然估计呢?下面通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法,由此可以清楚地看出EM算法的作用。
我们面对一个含有隐变量的概率模型,目标是极大化观测数据(不完全数据)Y关于参数 的对数似然函数,即极大化:
L
(
θ
)
=
l
o
g
P
(
Y
∣
θ
)
=
l
o
g
∑
z
P
(
Y
,
Z
∣
θ
)
=
l
o
g
(
∑
z
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
)
L(\theta)=log P(Y|\theta)=log \sum_z P(Y,Z|\theta)\\=log(\sum_zP(Y|Z,\theta)P(Z|\theta))
L(θ)=logP(Y∣θ)=logz∑P(Y,Z∣θ)=log(z∑P(Y∣Z,θ)P(Z∣θ))
注意到这一极大化的主要困难是上式中有未观测数据并有包含和(或积分)的对数。
事实上,EM算法是通过迭代逐步近似极大化
L
(
θ
)
L(\theta)
L(θ)的。假设在第i次迭代后θ的估计值θ(i)。我们希望新估计值θ能使L(θ)增加,即L(θ)>L(θ(i)),并逐步达到极大值。为此,考虑二者的差:
L
(
θ
)
−
L
(
θ
(
i
)
)
=
l
o
g
(
∑
z
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
)
−
l
o
g
P
(
Y
∣
θ
(
i
)
)
L(\theta)-L(\theta^{(i)})=log(\sum_zP(Y|Z,\theta)P(Z|\theta))-log P(Y|\theta^{(i)})
L(θ)−L(θ(i))=log(z∑P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θ(i))
利用Jensen不等式的得到其下界:
L ( θ ) − L ( θ ( i ) ) = log ( ∑ Z P ( Y ∣ Z , θ ( i ) ) P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Y ∣ Z , θ ( i ) ) ) − log P ( Y ∣ θ ( i ) ) ⩾ ∑ Z P ( Z ∣ Y , θ ( i ) ) log P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) − log P ( Y ∣ θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) log P ( Y ∣ Z , θ ) P ( Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) P ( Y ∣ θ ( i ) ) \begin{aligned} L(\theta)-L\left(\theta^{(i)}\right) &=\log \left(\sum_{Z} P\left(Y \mid Z, \theta^{(i)}\right) \frac{P(Y \mid Z, \theta) P(Z \mid \theta)}{P\left(Y \mid Z, \theta^{(i)}\right)}\right)-\log P\left(Y \mid \theta^{(i)}\right) \\ & \geqslant \sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log \frac{P(Y \mid Z, \theta) P(Z \mid \theta)}{P\left(Z \mid Y, \theta^{(i)}\right)}-\log P\left(Y \mid \theta^{(i)}\right) \\ &=\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log \frac{P(Y \mid Z, \theta) P(Z \mid \theta)}{P\left(Z \mid Y, \theta^{(i)}\right) P\left(Y \mid \theta^{(i)}\right)} \\ \end{aligned} L(θ)−L(θ(i))=log(Z∑P(Y∣Z,θ(i))P(Y∣Z,θ(i))P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θ(i))⩾Z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣Z,θ)P(Z∣θ)−logP(Y∣θ(i))=Z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ)
令
B
(
θ
,
θ
(
i
)
)
=
^
L
(
θ
(
i
)
)
+
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
P
(
Y
∣
θ
(
i
)
)
B\left(\theta, \theta^{(i)}\right) \hat{=} L\left(\theta^{(i)}\right)+\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log \frac{P(Y \mid Z, \theta) P(Z \mid \theta)}{P\left(Z \mid Y, \theta^{(i)}\right) P\left(Y \mid \theta^{(i)}\right)}
B(θ,θ(i))=^L(θ(i))+Z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ)
记上式为(2)式
则
L
(
θ
)
≥
B
(
θ
,
θ
(
i
)
)
L(\theta) \ge B(\theta,\theta^{(i)})
L(θ)≥B(θ,θ(i))
即函数
B
(
θ
,
θ
(
i
)
)
B(\theta,\theta^{(i)})
B(θ,θ(i))是
L
(
θ
)
L(\theta)
L(θ)的一个下界,而且由(2)式可知,
L
(
θ
(
i
)
)
=
B
(
θ
(
i
)
,
θ
(
i
)
)
L(\theta^{(i)}) = B(\theta^{(i)},\theta^{(i)})
L(θ(i))=B(θ(i),θ(i))
因此,任何可以使
B
(
θ
,
θ
(
i
)
)
B(\theta,\theta^{(i)})
B(θ,θ(i))增大的
θ
\theta
θ,也可以使L(θ)增大。为了使L(θ)有尽可能大的增大,选择θ(i+1)使得
B
(
θ
,
θ
(
i
)
)
B(\theta,\theta^{(i)})
B(θ,θ(i))达到极大,即
θ
i
+
1
=
a
r
g
m
a
x
B
(
θ
,
θ
(
i
)
)
\theta^{i+1} = argmaxB(\theta,\theta^{(i)})
θi+1=argmaxB(θ,θ(i))
记上述式子为(3)式
现在求θ(i+1)的表达式,省去对θ的极大化而言是常数的项,由(1)、(2)、(3)式可知:
θ
(
i
+
1
)
=
arg
max
θ
(
L
(
θ
(
i
)
)
+
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
P
(
Y
∣
θ
(
i
)
)
)
=
arg
max
θ
(
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
(
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
)
)
=
arg
max
θ
(
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
,
Z
∣
θ
)
)
=
arg
max
θ
Q
(
θ
,
θ
(
i
)
)
\begin{aligned} \theta^{(i+1)} &=\arg \max _{\theta}\left(L\left(\theta^{(i)}\right)+\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log \frac{P(Y \mid Z, \theta) P(Z \mid \theta)}{P\left(Z \mid Y, \theta^{(i)}\right) P\left(Y \mid \theta^{(i)}\right)}\right) \\ &=\arg \max _{\theta}\left(\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log (P(Y \mid Z, \theta) P(Z \mid \theta))\right) \\ &=\arg \max _{\theta}\left(\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log P(Y, Z \mid \theta)\right) \\ &=\arg \max _{\theta} Q\left(\theta, \theta^{(i)}\right) \end{aligned}
θ(i+1)=argθmax(L(θ(i))+Z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ))=argθmax(Z∑P(Z∣Y,θ(i))log(P(Y∣Z,θ)P(Z∣θ)))=argθmax(Z∑P(Z∣Y,θ(i))logP(Y,Z∣θ))=argθmaxQ(θ,θ(i))
上述式子等价于EM算法的一次迭代,即求Q函数及其极大化,EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
上图给出了EM算法的直观解释,图中上方曲线为 L ( θ ) L(\theta) L(θ),下方曲线为 B ( θ , θ ( i ) ) B(\theta,\theta^{(i)}) B(θ,θ(i))。由前面推导可以知道, B ( θ , θ ( i ) ) B(\theta,\theta^{(i)}) B(θ,θ(i))为对数似然函数 L ( θ ) L(\theta) L(θ)的下界。
两个函数在点 θ = θ ( i ) \theta=\theta^{(i)} θ=θ(i)处相等。
EM算法找到下一个点 θ ( i ) \theta^{(i)} θ(i)使函数 B ( θ , θ ( i ) ) B(\theta,\theta^{(i)}) B(θ,θ(i))极大化,也使函数 Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q(θ,θ(i))极大化。这时由于 L ( θ ) ≥ B ( θ , θ ( i ) ) L(\theta) \ge B(\theta,\theta^{(i)}) L(θ)≥B(θ,θ(i)),函数 B ( θ , θ ( i ) ) B(\theta,\theta^{(i)}) B(θ,θ(i))的增加,保证对数似然函数 L ( θ ) L(\theta) L(θ)在每次迭代中也是增加的。EM算法在点 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)重新计算Q函数 值,进行下一次迭代。在这个过程中,对数似然函数L(θ)不断增大。 从图可以推断出EM算法不能保证找到全局最优值。