EM算法的本质:
EM算法其实也是利用了极大似然估计的思想,只不过是极大似然估计是在观测样本不缺失的情况下,找到模型参数使得观测样本出现的概率最大。而EM算法要解决的问题是,当观测样本数据不完整(缺失)的情况下,找到模型参数使得观测样本出现的概率最大并且补全缺失的数据。与传统的极大似然估计求参数的方法不同,此时因为含有未知的数据,就不能直接极大化对数似然函数来求解模型参数。这时候EM算法就派上用场了。
EM算法概述:
EM算法解决这个问题使用的是迭代的方法,既然我们无法直接求解出模型的参数,那就先猜想缺失的数据(EM算法的E步),然后再利用观测样本和猜想的缺失数据极大化对数似然函数,求解模型的参数(EM的M步)。由于我们之前的缺失的数据是猜想来的,故第一次求解出来的参数大概率不是最终的结果,不过没关系,基于当前的参数,我们可以继续利用样本猜测缺失的数据(EM算法的E步),然后继续极大化对数似然,求解我们的模型参数(EM算法的M步)。以此类推,不断的迭代下去,直到模型分布参数基本无变化,算法收敛,找到合适的模型参数。
K-Means与EM算法的联系:
一个最直观的例子就是,之前提到过的K-Means聚类算法,EM算法和K-Means聚类算法都是迭代求解的过程,在K-Means中,每个样本所属的类就可以看成是一个隐变量,k个中心点即为我们要求解的参数。首先要初始化k个中心点,然后使得总体类内距离最小,即把样本指派到与其最近的中心所属的类中,得到聚类(E step),然后根据得到的聚类重新确定k个中心点(M step)。K-Means聚类算法
EM算法的推导
对于m个样本的观测数据
x
=
(
x
(
1
)
,
x
(
2
)
,
.
.
.
x
(
m
)
)
x=(x^{(1)},x^{(2)},...x^{(m)})
x=(x(1),x(2),...x(m)),极大化模型的对数似然函数,找到对于的参数
θ
=
a
r
g
max
θ
∑
i
=
1
m
l
o
g
P
(
x
(
i
)
;
θ
)
\theta = arg \max \limits_{\theta}\sum\limits_{i=1}^m logP(x^{(i)};\theta)
θ=argθmaxi=1∑mlogP(x(i);θ)
但是,如果我们有未知的样本
z
=
(
z
(
1
)
,
z
(
2
)
,
.
.
.
z
(
m
)
)
z=(z^{(1)},z^{(2)},...z^{(m)})
z=(z(1),z(2),...z(m)),就要重新考虑我们的对数似然函数:
θ
=
a
r
g
max
θ
∑
i
=
1
m
l
o
g
P
(
x
(
i
)
;
θ
)
=
a
r
g
max
θ
∑
i
=
1
m
l
o
g
∑
z
(
i
)
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
\theta = arg \max \limits_{\theta}\sum\limits_{i=1}^m logP(x^{(i)};\theta) = arg \max \limits_{\theta}\sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta)
θ=argθmaxi=1∑mlogP(x(i);θ)=argθmaxi=1∑mlogz(i)∑P(x(i),z(i);θ)
此时,我们可以利用Jensen不等式变形为:
∑
i
=
1
m
l
o
g
∑
z
(
i
)
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
=
∑
i
=
1
m
l
o
g
∑
z
(
i
)
Q
i
(
z
(
i
)
)
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
Q
i
(
z
(
i
)
)
≥
∑
i
=
1
m
∑
z
(
i
)
Q
i
(
z
(
i
)
)
l
o
g
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
Q
i
(
z
(
i
)
)
(1)
\begin{align} \sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta) & = \sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)}, z^{(i)};\theta)}{Q_i(z^{(i)})} \\ & \geq \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)}, z^{(i)};\theta)}{Q_i(z^{(i)})} \end{align}\tag{1}
i=1∑mlogz(i)∑P(x(i),z(i);θ)=i=1∑mlogz(i)∑Qi(z(i))Qi(z(i))P(x(i),z(i);θ)≥i=1∑mz(i)∑Qi(z(i))logQi(z(i))P(x(i),z(i);θ)(1)
其中
Q
i
(
z
(
i
)
)
Q_i(z^{(i)})
Qi(z(i))代表丢失变量z的概率。
Jensen不等式是这样定义的:
f
(
E
(
x
)
)
≥
E
(
f
(
x
)
)
如果
f
(
x
)
是凹函数
(2)
f(E(x)) \geq E(f(x))\;\; 如果f(x) 是凹函数 \tag{2}
f(E(x))≥E(f(x))如果f(x)是凹函数(2)
上面公式的f(x)就是log函数。
此时如果要满足Jensen不等式的等号,当且仅当:(2)式中的x为常数,对应到公式(1)中,即:
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
Q
i
(
z
(
i
)
)
=
c
,
c
为常数
\frac{P(x^{(i)}, z^{(i)};\theta)}{Q_i(z^{(i)})} =c, c为常数
Qi(z(i))P(x(i),z(i);θ)=c,c为常数
因为:
∑
z
Q
i
(
z
(
i
)
)
=
1
\sum\limits_{z}Q_i(z^{(i)}) =1
z∑Qi(z(i))=1
从上面两式,我们可以得到:
Q
i
(
z
(
i
)
)
=
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
∑
z
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
=
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
P
(
x
(
i
)
;
θ
)
=
P
(
z
(
i
)
∣
x
(
i
)
;
θ
)
)
Q_i(z^{(i)}) = \frac{P(x^{(i)}, z^{(i)};\theta)}{\sum\limits_{z}P(x^{(i)}, z^{(i)};\theta)} = \frac{P(x^{(i)}, z^{(i)};\theta)}{P(x^{(i)};\theta)} = P( z^{(i)}|x^{(i)};\theta))
Qi(z(i))=z∑P(x(i),z(i);θ)P(x(i),z(i);θ)=P(x(i);θ)P(x(i),z(i);θ)=P(z(i)∣x(i);θ))
到此,我们便求出了丢失变量
z
z
z的概率。
有了 Q i ( z ( i ) ) Q_i(z^{(i)}) Qi(z(i)),怎么更新模型的参数使得对数似然函数逐渐增加呢?(M-step)
我们需要最大化似然函数,即(1)式,也就是最大化(1)式中的下界,把其中的常数想去掉,即最大化:
a
r
g
max
θ
∑
i
=
1
m
∑
z
(
i
)
Q
i
(
z
(
i
)
)
l
o
g
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
arg \max \limits_{\theta} \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}, z^{(i)};\theta)}
argθmaxi=1∑mz(i)∑Qi(z(i))logP(x(i),z(i);θ)
知道了
M
−
s
t
e
p
M-step
M−step,那么
E
−
s
t
e
p
E-step
E−step呢,
E
−
s
t
e
p
E-step
E−step即计算联合分布的条件概率期望:
Q
i
(
z
(
i
)
)
=
P
(
z
(
i
)
∣
x
(
i
)
,
θ
j
)
)
Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)},\theta^{j}))
Qi(z(i))=P(z(i)∣x(i),θj))
L ( θ , θ j ) = ∑ i = 1 m ∑ z ( i ) Q i ( z ( i ) ) l o g P ( x ( i ) , z ( i ) ; θ ) (3) L(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}, z^{(i)};\theta)}\tag{3} L(θ,θj)=i=1∑mz(i)∑Qi(z(i))logP(x(i),z(i);θ)(3)
EM算法的流程:
输入观测数据 x = ( x ( 1 ) , x ( 2 ) , . . . x ( m ) ) x=(x^{(1)},x^{(2)},...x^{(m)}) x=(x(1),x(2),...x(m)), z = ( z ( 1 ) , z ( 2 ) , . . . z ( m ) ) z=(z^{(1)},z^{(2)},...z^{(m)}) z=(z(1),z(2),...z(m)),联合分布 p ( x , z ; θ ) p(x,z ;\theta) p(x,z;θ),条件分布 p ( z ∣ x ; θ ) p(z|x; \theta) p(z∣x;θ),最大迭代次数 J J J,
-
初始化模型的参数 θ \theta θ为 θ 0 \theta_0 θ0
-
for j from 1 to J开始EM算法迭代:
E -step:计算联合分布的条件概率期望
Q i ( z ( i ) ) = P ( z ( i ) ∣ x ( i ) , θ j ) ) L ( θ , θ j ) = ∑ i = 1 m ∑ z ( i ) Q i ( z ( i ) ) l o g P ( x ( i ) , z ( i ) ; θ ) Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)},\theta^{j}))\\ L(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}, z^{(i)};\theta)} Qi(z(i))=P(z(i)∣x(i),θj))L(θ,θj)=i=1∑mz(i)∑Qi(z(i))logP(x(i),z(i);θ)
M-step:极大化
L
(
θ
,
θ
j
)
L(\theta, \theta^{j})
L(θ,θj),得到
θ
j
+
1
\theta^{j+1}
θj+1。
θ
j
+
1
=
a
r
g
max
θ
L
(
θ
,
θ
j
)
\theta^{j+1} = arg \max \limits_{\theta}L(\theta, \theta^{j})
θj+1=argθmaxL(θ,θj)
如果
θ
j
+
1
\theta^{j+1}
θj+1 收敛,则停止迭代,否则回到E-step继续迭代。
EM算法收敛性的思考:
(1)EM算法能够收敛吗?
(2)EM算法能够收敛到全局最大值么?
第一个问题:EM算法要想要收敛而且想要达到最大值,则必须 保证EM算法迭代过程中对数似然函数是一直增大即:
∑
i
=
1
m
l
o
g
P
(
x
(
i
)
;
θ
j
+
1
)
≥
∑
i
=
1
m
l
o
g
P
(
x
(
i
)
;
θ
j
)
\sum\limits_{i=1}^m logP(x^{(i)};\theta^{j+1}) \geq \sum\limits_{i=1}^m logP(x^{(i)};\theta^{j})
i=1∑mlogP(x(i);θj+1)≥i=1∑mlogP(x(i);θj)
由于:
L
(
θ
,
θ
j
)
=
∑
i
=
1
m
∑
z
(
i
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
)
)
l
o
g
P
(
x
(
i
)
,
z
(
i
)
;
θ
)
L(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j}))log{P(x^{(i)}, z^{(i)};\theta)}
L(θ,θj)=i=1∑mz(i)∑P(z(i)∣x(i);θj))logP(x(i),z(i);θ)
令:
H
(
θ
,
θ
j
)
=
∑
i
=
1
m
∑
z
(
i
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
)
)
l
o
g
P
(
z
(
i
)
∣
x
(
i
)
;
θ
)
H(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j}))log{P( z^{(i)}|x^{(i)};\theta)}
H(θ,θj)=i=1∑mz(i)∑P(z(i)∣x(i);θj))logP(z(i)∣x(i);θ)
上面两式子相减得到:
∑
i
=
1
m
l
o
g
P
(
x
(
i
)
;
θ
)
=
L
(
θ
,
θ
j
)
−
H
(
θ
,
θ
j
)
\sum\limits_{i=1}^m logP(x^{(i)};\theta) = L(\theta, \theta^{j}) - H(\theta, \theta^{j})
i=1∑mlogP(x(i);θ)=L(θ,θj)−H(θ,θj)
在上式中分别取θ为
θ
j
θ^j
θj和
θ
j
+
1
θ^{j+1}
θj+1,并相减得到:
∑
i
=
1
m
l
o
g
P
(
x
(
i
)
;
θ
j
+
1
)
−
∑
i
=
1
m
l
o
g
P
(
x
(
i
)
;
θ
j
)
=
[
L
(
θ
j
+
1
,
θ
j
)
−
L
(
θ
j
,
θ
j
)
]
−
[
H
(
θ
j
+
1
,
θ
j
)
−
H
(
θ
j
,
θ
j
)
]
\sum\limits_{i=1}^m logP(x^{(i)};\theta^{j+1}) - \sum\limits_{i=1}^m logP(x^{(i)};\theta^{j}) = [L(\theta^{j+1}, \theta^{j}) - L(\theta^{j}, \theta^{j}) ] -[H(\theta^{j+1}, \theta^{j}) - H(\theta^{j}, \theta^{j}) ]
i=1∑mlogP(x(i);θj+1)−i=1∑mlogP(x(i);θj)=[L(θj+1,θj)−L(θj,θj)]−[H(θj+1,θj)−H(θj,θj)]
要证明EM算法的收敛性,我们只需要证明上式的右边是非负的即可。
由于
θ
j
+
1
\theta^{j+1}
θj+1使得
L
(
θ
,
θ
j
)
L(\theta, \theta^{j})
L(θ,θj)极大,所以
L
(
θ
j
+
1
,
θ
j
)
−
L
(
θ
j
,
θ
j
)
≥
0
L(\theta^{j+1}, \theta^{j}) - L(\theta^{j}, \theta^{j}) \geq 0
L(θj+1,θj)−L(θj,θj)≥0
所以只需要证明右边式子小于0即可:
H
(
θ
j
+
1
,
θ
j
)
−
H
(
θ
j
,
θ
j
)
=
∑
i
=
1
m
∑
z
(
i
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
)
l
o
g
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
+
1
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
)
≤
∑
i
=
1
m
l
o
g
(
∑
z
(
i
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
+
1
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
)
)
=
∑
i
=
1
m
l
o
g
(
∑
z
(
i
)
P
(
z
(
i
)
∣
x
(
i
)
;
θ
j
+
1
)
)
=
0
\begin{align} H(\theta^{j+1}, \theta^{j}) - H(\theta^{j}, \theta^{j}) & = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j})log\frac{P( z^{(i)}|x^{(i)};\theta^{j+1})}{P( z^{(i)}|x^{(i)};\theta^j)} \\ & \leq \sum\limits_{i=1}^mlog(\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j})\frac{P( z^{(i)}|x^{(i)};\theta^{j+1})}{P( z^{(i)}|x^{(i)};\theta^j)}) \\ & = \sum\limits_{i=1}^mlog(\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j+1})) = 0 \end{align}
H(θj+1,θj)−H(θj,θj)=i=1∑mz(i)∑P(z(i)∣x(i);θj)logP(z(i)∣x(i);θj)P(z(i)∣x(i);θj+1)≤i=1∑mlog(z(i)∑P(z(i)∣x(i);θj)P(z(i)∣x(i);θj)P(z(i)∣x(i);θj+1))=i=1∑mlog(z(i)∑P(z(i)∣x(i);θj+1))=0
其中第(4)式用到了Jensen不等式,只不过和第二节的使用相反而已,第(5)式用到了概率分布累积为1的性质。
至此,我们得到了: ∑ i = 1 m l o g P ( x ( i ) ; θ j + 1 ) − ∑ i = 1 m l o g P ( x ( i ) ; θ j ) ≥ 0 \sum\limits_{i=1}^m logP(x^{(i)};\theta^{j+1}) - \sum\limits_{i=1}^m logP(x^{(i)};\theta^{j}) \geq 0 i=1∑mlogP(x(i);θj+1)−i=1∑mlogP(x(i);θj)≥0。 证明了EM算法的收敛性。
从上面的推导可以看出,EM算法可以保证收敛到一个稳定点,但是却不能保证收敛到全局的极大值点,因此它是局部最优的算法,当然,如果我们的优化目标 L ( θ , θ j ) L(\theta, \theta^{j}) L(θ,θj)是凸的,则EM算法可以保证收敛到全局最大值,这点和梯度下降法这样的迭代算法相同。