期望极大算法(EM算法)
看本篇之前如果不了解似然函数和极大似然估计的话可能有点要命,请先弄清楚相应的知识。另外,如果有符号/名词的含义不太明白,可拉到文章的最下方看看相应的说明。
在极大似然估计的博客中,我们在公式层面上讨论了怎么进行概率模型参数的极大似然估计。在文章中举的掷色子的例子中(不含有隐变量),我们也了解到可以通过求解函数的导数求得极值来得到这个极大似然估计。
在本篇文章中,我们将探讨含有隐变量的概率模型参数的极大似然估计算法—EM算法。
隐变量
首先要明白啥是隐变量,看到李航老师《统计学习方法》中举的例子:
还是抛硬币:假设有3枚硬币,分别记作A,B,C。这些硬币正面出现的概率分别是 π \pi π, p p p和 q q q.进行如下掷硬币试验:先掷硬币A,根据其结果选出硬币B或硬币C,正面选硬币B,反面选硬币C;然后掷选出的硬币,掷硬币的结果,出现正面记作1,出现反面记作0;独立地重复n次试验(n=10),观测结果如下:1,1,0,1,0,0,1,0,1,1
假设只能观测到掷硬币的结果,不能观测掷硬币的过程。问如何估计三硬币正面出现的概率,即三硬币模型的参数。
在这个例子中,我们不知道第一次抛硬币的结果,于是第二次用的是哪个硬币抛的也不知道,而这个不知道的变量又影响着最终的观察,因此我们将这种影响结果而又无法观测到的变量称之为隐变量。
带隐变量的似然函数
由于带有隐变量,我们需要将其可能造成的影响全部考虑进来,因此硬币模型就得写成这样:
P
(
y
∣
θ
)
=
∑
z
P
(
y
,
z
∣
θ
)
=
∑
z
P
(
z
∣
θ
)
P
(
y
∣
z
,
θ
)
=
π
p
y
(
1
−
p
)
1
−
y
+
(
1
−
π
)
q
y
(
1
−
q
)
1
−
y
\begin{aligned} P(y \mid \theta) &= \sum_{z} P(y, z \mid \theta)=\sum_{z} P(z \mid \theta) P(y \mid z, \theta) \\& =\pi p^{y}(1-p)^{1-y}+(1-\pi) q^{y}(1-q)^{1-y} \end{aligned}
P(y∣θ)=z∑P(y,z∣θ)=z∑P(z∣θ)P(y∣z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y
这里,随机变量
y
y
y是观测变量,表示一次试验观测的结果是1或0;随机变量
z
z
z是隐变量,表示未观测到的掷硬币А的结果:
θ
=
(
π
,
p
,
q
)
\theta=(\pi,p,q)
θ=(π,p,q)是模型参数,这一模型是以上数据的生成模型(随机变量
y
y
y的数据可以观测,随机变量
z
z
z的数据不可观测)。
类比抛掷一个硬币的模型:
P
(
y
∣
θ
)
=
θ
y
(
1
−
θ
)
1
−
y
\begin{aligned} P(y \mid \theta) & = \theta^{y}(1-\theta)^{1-y} \end{aligned}
P(y∣θ)=θy(1−θ)1−y
显然上面的模型要多一个隐含项。
有了硬币的模型,将观测数据表示为
Y
=
(
Y
1
,
Y
2
,
⋯
,
Y
n
)
T
Y=\left(Y_{1}, Y_{2}, \cdots, Y_{n}\right)^{\mathrm{T}}
Y=(Y1,Y2,⋯,Yn)T,未观测数据表示为
Z
=
(
Z
1
,
Z
2
,
⋯
,
Z
n
)
T
Z=\left(Z_{1}, Z_{2}, \cdots, Z_{n}\right)^{\mathrm{T}}
Z=(Z1,Z2,⋯,Zn)T,则李航老师的书告诉我们观测数据的似然函数为:
P
(
Y
∣
θ
)
=
∑
Z
P
(
Z
∣
θ
)
P
(
Y
∣
Z
,
θ
)
P(Y \mid \theta)=\sum_{Z} P(Z \mid \theta) P(Y \mid Z, \theta)
P(Y∣θ)=Z∑P(Z∣θ)P(Y∣Z,θ)
这里本人看了半天才看明白。因为在极大似然估计的博客中我们知道各观测值之间应该是一个连乘的关系,怎么这里弄出的是个连加的符号,而连乘符号在下一道式子中突然又出现了?
这里我们将式子写成下面这样就好懂很多了:
P
(
Y
∣
θ
)
=
∏
j
=
1
n
∑
z
j
P
(
z
j
∣
θ
)
P
(
y
j
∣
z
j
,
θ
)
=
∑
Z
P
(
Z
∣
θ
)
P
(
Y
∣
Z
,
θ
)
P(Y \mid \theta)=\prod_{j=1}^{n}\sum_{z_j} P(z_j \mid \theta) P(y_j \mid z_j, \theta) =\sum_{Z} P(Z \mid \theta) P(Y \mid Z, \theta)
P(Y∣θ)=j=1∏nzj∑P(zj∣θ)P(yj∣zj,θ)=Z∑P(Z∣θ)P(Y∣Z,θ)
它确实是连乘了的,只是李老师的书把序列当成一个整体写了,
P
(
Z
∣
θ
)
P(Z \mid \theta)
P(Z∣θ)表示的是在
θ
\theta
θ的模型参数下,序列
Z
Z
Z整体出现的概率,
P
(
Y
∣
Z
,
θ
)
P(Y \mid Z, \theta)
P(Y∣Z,θ)也同理,连乘符号隐含在其中。
等式最终可以写成:
P
(
Y
∣
θ
)
=
∏
j
=
1
n
[
π
p
y
j
(
1
−
p
)
1
−
y
j
+
(
1
−
π
)
q
y
j
(
1
−
q
)
1
−
y
j
]
P(Y \mid \theta)=\prod_{j=1}^{n}\left[\pi p^{y_{j}}(1-p)^{1-y_{j}}+(1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}\right]
P(Y∣θ)=j=1∏n[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj]
即得到了抛硬币的似然函数。
带隐变量的极大似然估计
为了求得硬币的具体参数
θ
=
(
π
,
p
,
q
)
\theta=(\pi, p, q)
θ=(π,p,q),我们需要做极大似然估计,结合极大似然估计的推导,我们知道实际上我们是要去解如下公式:
θ
^
=
arg
max
θ
log
P
(
Y
∣
θ
)
\hat{\theta}=\arg \max _{\theta} \log P(Y \mid \theta)
θ^=argθmaxlogP(Y∣θ)
但这里就出现了最要命的问题,由于公式中变量太多,没有解析解!
EM算法导出
为了求解上面的函数,EM算法就出现了!在前面的推导中我们知道,要求带隐变量的极大似然估计,实际上是要极大化
L
(
θ
)
=
log
P
(
Y
∣
θ
)
=
log
(
∑
Z
P
(
Y
,
Z
∣
θ
)
)
=
log
(
∑
Z
P
(
Z
∣
θ
)
P
(
Y
∣
Z
,
θ
)
)
L(\theta)=\log P(Y \mid \theta)=\log\left(\sum_{Z} P(Y, Z \mid \theta) \right)\\=\log\left(\sum_{Z} P(Z \mid \theta) P(Y \mid Z, \theta)\right)
L(θ)=logP(Y∣θ)=log(Z∑P(Y,Z∣θ))=log(Z∑P(Z∣θ)P(Y∣Z,θ))
而这一极大化因为包含了未观测数据且有和的对数等等,无法直接求得解析解。
为了解决这个问题,EM算法提出通过迭代逐步近似极大化 L ( θ ) L(\theta) L(θ)的操作。操作的具体逻辑如下:
假设在第
i
i
i次迭代后
θ
\theta
θ的估计值是
θ
(
i
)
\theta^{(i)}
θ(i),那么出于极大化的考虑,我们希望新的估计值
θ
\theta
θ要使
L
(
θ
)
L(\theta)
L(θ)增加,即
L
(
θ
)
>
L
(
θ
(
i
)
)
L(\theta)>L(\theta^{(i)})
L(θ)>L(θ(i)),这样反复迭代最后达到极大值。为此,考虑两者的差:
L
(
θ
)
−
L
(
θ
(
i
)
)
=
log
(
∑
Z
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
)
−
log
P
(
Y
∣
θ
(
i
)
)
L(\theta)-L\left(\theta^{(i)}\right)=\log \left(\sum_{Z} P(Y \mid Z, \theta) P(Z \mid \theta)\right)-\log P\left(Y \mid \theta^{(i)}\right)
L(θ)−L(θ(i))=log(Z∑P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θ(i))
利用琴生(jensen)不等式加权形式:
f
(
∑
j
λ
j
y
j
)
⩾
∑
j
λ
j
f
(
y
j
)
,
其
中
λ
j
⩾
0
,
∑
j
λ
j
=
1
f(\sum_{j}\lambda_{j}y_{j})\geqslant\sum_{j}\lambda_{j}f(y_{j}),其中\lambda_{j}\geqslant0,\sum_{j}\lambda_{j}=1
f(j∑λjyj)⩾j∑λjf(yj),其中λj⩾0,j∑λj=1
取
f
(
y
j
)
=
l
o
g
y
j
f(y_j)=logy_j
f(yj)=logyj,有
log
∑
j
λ
j
y
j
⩾
∑
j
λ
j
log
y
j
,
其
中
λ
j
⩾
0
,
∑
j
λ
j
=
1
\log\sum_{j}\lambda_{j}y_{j}\geqslant\sum_{j}\lambda_{j}\log y_{j},其中\lambda_{j}\geqslant0,\sum_{j}\lambda_{j}=1
logj∑λjyj⩾j∑λjlogyj,其中λj⩾0,j∑λj=1
则:
L
(
θ
)
−
L
(
θ
(
i
)
)
=
log
(
∑
Z
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
)
−
log
P
(
Y
∣
θ
(
i
)
)
=
log
(
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
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
)
)
−
log
P
(
Y
∣
θ
(
i
)
)
L(\theta)-L\left(\theta^{(i)}\right)=\log \left(\sum_{Z} P(Y \mid Z, \theta) P(Z \mid \theta)\right)-\log P\left(Y \mid \theta^{(i)}\right)\\ =\log \left(\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \frac{P(Y \mid Z, \theta) P(Z \mid \theta)}{P\left(Z \mid Y, \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)
L(θ)−L(θ(i))=log(Z∑P(Y∣Z,θ)P(Z∣θ))−logP(Y∣θ(i))=log(Z∑P(Z∣Y,θ(i))P(Z∣Y,θ(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
)
)
=
1
\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right)=1
∑ZP(Z∣Y,θ(i))=1且
log
P
(
Y
∣
θ
(
i
)
)
\log P\left(Y \mid \theta^{(i)}\right)
logP(Y∣θ(i))不存在
Z
Z
Z变量,因此
log
P
(
Y
∣
θ
(
i
)
)
=
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
∣
θ
(
i
)
)
\log P\left(Y \mid \theta^{(i)}\right)=\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right)\log P\left(Y \mid \theta^{(i)}\right)
logP(Y∣θ(i))=Z∑P(Z∣Y,θ(i))logP(Y∣θ(i))
有:
L
(
θ
)
−
L
(
θ
(
i
)
)
⩾
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
P
(
Y
∣
θ
(
i
)
)
L(\theta)-L\left(\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)P\left(Y \mid \theta^{(i)}\right)}
L(θ)−L(θ(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∣θ)
得到
L
(
θ
)
⩾
B
(
θ
,
θ
(
i
)
)
L(\theta)\geqslant B\left(\theta, \theta^{(i)}\right)
L(θ)⩾B(θ,θ(i))
即函数
B
(
θ
,
θ
(
i
)
)
B\left(\theta, \theta^{(i)}\right)
B(θ,θ(i))是
L
(
θ
)
L(\theta)
L(θ)的一个下界,当取
θ
=
θ
(
i
)
\theta=\theta^{(i)}
θ=θ(i),有
B
(
θ
(
i
)
,
θ
(
i
)
)
=
L
(
θ
(
i
)
)
+
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
∣
Z
,
θ
(
i
)
)
P
(
Z
∣
θ
(
i
)
)
P
(
Z
∣
Y
,
θ
(
i
)
)
P
(
Y
∣
θ
(
i
)
)
=
L
(
θ
(
i
)
)
+
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
1
=
L
(
θ
(
i
)
)
B\left(\theta^{(i)}, \theta^{(i)}\right)=L\left(\theta^{(i)}\right)+\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log \frac{P(Y \mid Z, \theta^{(i)}) P(Z \mid \theta^{(i)})}{P\left(Z \mid Y, \theta^{(i)}\right) P\left(Y \mid \theta^{(i)}\right)}\\=L\left(\theta^{(i)}\right)+\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log 1=L\left(\theta^{(i)}\right)
B(θ(i),θ(i))=L(θ(i))+Z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ(i))P(Z∣θ(i))=L(θ(i))+Z∑P(Z∣Y,θ(i))log1=L(θ(i))
也就是说
B
(
θ
,
θ
(
i
)
)
B\left(\theta, \theta^{(i)}\right)
B(θ,θ(i))至少是个极大值为
L
(
θ
(
i
)
)
L\left(\theta^{(i)}\right)
L(θ(i))的函数,只要每次迭代我们选取
θ
\theta
θ使得
B
(
θ
,
θ
(
i
)
)
B\left(\theta, \theta^{(i)}\right)
B(θ,θ(i))达到极大值,则我们至少得到的是一个不差于之前的似然估计,如果迭代没收敛,我们则将获得更好的似然估计。
现在问题转换为选择
θ
(
i
+
1
)
\theta^{(i+1)}
θ(i+1)使得
B
(
θ
,
θ
(
i
)
)
B\left(\theta, \theta^{(i)}\right)
B(θ,θ(i))获得极大值,即
θ
(
i
+
1
)
=
arg
max
θ
B
(
θ
,
θ
(
i
)
)
=
arg
max
θ
(
L
(
θ
(
i
)
)
+
∑
z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
P
(
Y
∣
θ
(
i
)
)
)
\theta^{(i+1)}=\arg \max _{\theta} B\left(\theta, \theta^{(i)}\right)\\=\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)
θ(i+1)=argθmaxB(θ,θ(i))=argθmax(L(θ(i))+z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ))
由于
L
(
θ
(
i
)
)
L\left(\theta^{(i)}\right)
L(θ(i))是定值,因此可以忽略,公式简化为:
arg
max
θ
(
∑
z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
∣
Z
,
θ
)
P
(
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
P
(
Y
∣
θ
(
i
)
)
)
\arg \max _{\theta}\left(\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(z∑P(Z∣Y,θ(i))logP(Z∣Y,θ(i))P(Y∣θ(i))P(Y∣Z,θ)P(Z∣θ))
又
P
(
Z
∣
Y
,
θ
(
i
)
)
P
(
Y
∣
θ
(
i
)
)
=
P
(
Z
,
Y
∣
θ
(
i
)
)
=
M
P\left(Z \mid Y, \theta^{(i)}\right) P\left(Y \mid \theta^{(i)}\right)=P\left(Z , Y\mid \theta^{(i)}\right)=M
P(Z∣Y,θ(i))P(Y∣θ(i))=P(Z,Y∣θ(i))=M,
M
M
M为常数,因此原式等于:
arg
max
θ
(
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
,
Z
∣
θ
)
)
=
arg
max
θ
Q
(
θ
,
θ
(
i
)
)
\arg \max _{\theta}\left(\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log P(Y, Z \mid \theta)\right) =\underset{\theta}{\arg \max } Q\left(\theta, \theta^{(i)}\right)
argθmax(Z∑P(Z∣Y,θ(i))logP(Y,Z∣θ))=θargmaxQ(θ,θ(i))
也就是说我们每次迭代就相当于要极大化一次
Q
(
θ
,
θ
(
i
)
)
Q\left(\theta, \theta^{(i)}\right)
Q(θ,θ(i))函数,经过足够次数的迭代,最终我们将得到一个相对准确的极大似然估计。
需要注意的是,由于EM是一个迭代算法,当迭代到 B ( θ , θ ( i ) ) B\left(\theta, \theta^{(i)}\right) B(θ,θ(i))的局部最大值点时,迭代可能会终止,因此不能保证最终结果收敛到极大值点。于是,初值的选取非常重要,常用的方法是选取几个不同的初值进行迭代,然后对得到的各个估计值加以比较,选择最好的一个。。。
EM算法实际操作步骤
上面我们已经完成了整个EM算法的推导,过程虽然坎坷,但最后的结论其实非常简单,那就是每次迭代求一次 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))的极大值。流程如下:
-
选择参数的初值 θ 0 \theta^{0} θ0,开始迭代;
-
E步:记 θ i \theta^{i} θi为第 i i i次迭代参数 θ \theta θ的估计值,在第 i + 1 i+1 i+1次迭代的E步,计算
Q ( θ , θ ( i ) ) = ∑ Z P ( Z ∣ Y , θ ( i ) ) log P ( Y , Z ∣ θ ) Q\left(\theta, \theta^{(i)}\right)=\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log P(Y, Z \mid \theta) Q(θ,θ(i))=Z∑P(Z∣Y,θ(i))logP(Y,Z∣θ)
这里, P ( Z ∣ Y , θ ( i ) ) P\left(Z \mid Y, \theta^{(i)}\right) P(Z∣Y,θ(i))是在给定观测数据 Y Y Y和当前的参数估计 θ ( i ) \theta^{(i)} θ(i)下隐变量数据 Z Z Z的条件概率分布; -
M步:求使 Q ( θ , θ ( i ) ) Q\left(\theta, \theta^{(i)}\right) Q(θ,θ(i))极大化的 θ \theta θ,确定第 i + 1 i+1 i+1次迭代的参数的估计值 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)。
-
重复2、3步直到收敛。
EM算法解决抛硬币问题
回到最初的抛硬币问题:
假设有3枚硬币,分别记作A,B,C。这些硬币正面出现的概率分别是 π \pi π, p p p和 q q q.进行如下掷硬币试验:先掷硬币A,根据其结果选出硬币B或硬币C,正面选硬币B,反面选硬币C;然后掷选出的硬币,掷硬币的结果,出现正面记作1,出现反面记作0;独立地重复n次试验(n=10),观测结果如下:1,1,0,1,0,0,1,0,1,1
假设只能观测到掷硬币的结果,不能观测掷硬币的过程。问如何估计三硬币正面出现的概率,即三硬币模型的参数。
在前面的推导中我们知道,要进行EM算法,首先要求得
Q
(
θ
,
θ
(
i
)
)
Q\left(\theta, \theta^{(i)}\right)
Q(θ,θ(i))函数,在这一具体问题中,有
Q
(
θ
,
θ
(
i
)
)
=
∑
Z
P
(
Z
∣
Y
,
θ
(
i
)
)
log
P
(
Y
,
Z
∣
θ
)
=
∑
Z
{
∏
j
=
1
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
log
[
∏
j
=
1
n
P
(
y
j
,
z
j
∣
θ
)
]
}
=
∑
Z
{
∏
j
=
1
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
∑
j
=
1
n
log
[
P
(
y
j
,
z
j
∣
θ
)
]
}
Q\left(\theta, \theta^{(i)}\right)=\sum_{Z} P\left(Z \mid Y, \theta^{(i)}\right) \log P(Y, Z \mid \theta)\\ =\sum_{Z}\left\{\prod_{j=1}^{n} P\left(z_j|y_j,\theta^{(i)}\right) \log\left[\prod_{j=1}^{n}P\left(y_j , z_j\mid \theta\right)\right] \right\}\\ =\sum_{Z}\left\{\prod_{j=1}^{n} P\left(z_j|y_j,\theta^{(i)}\right) \sum_{j=1}^{n}\log\left[P\left(y_j , z_j\mid \theta\right)\right] \right\}
Q(θ,θ(i))=Z∑P(Z∣Y,θ(i))logP(Y,Z∣θ)=Z∑{j=1∏nP(zj∣yj,θ(i))log[j=1∏nP(yj,zj∣θ)]}=Z∑{j=1∏nP(zj∣yj,θ(i))j=1∑nlog[P(yj,zj∣θ)]}
最开始本人以为这样大概就KO了,结果发现上式中有个非常让人崩溃的连乘符号,需要想办法把它弄掉,否则根本无法计算!这里采用这篇文章的方法,对$\prod_{j=1}^{n} P\left(z_j|y_j,\theta^{(i)}\right)
和
和
和\sum_{j=1}^{n}\log\left[P\left(y_j , z_j\mid \theta\right)\right]$进行分拆,先拆出第一次投币的结构,有
∏
j
=
1
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
=
P
(
z
1
∣
y
1
,
θ
(
i
)
)
×
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
\prod_{j=1}^{n} P\left(z_j|y_j,\theta^{(i)}\right) =P\left(z_1|y_1,\theta^{(i)}\right) × \prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)
j=1∏nP(zj∣yj,θ(i))=P(z1∣y1,θ(i))×j=2∏nP(zj∣yj,θ(i))
∑ j = 1 n log [ P ( y j , z j ∣ θ ) ] = log [ P ( y 1 , z 1 ∣ θ ) ] + ∑ j = 2 n log [ P ( y j , z j ∣ θ ) ] \sum_{j=1}^{n}\log\left[P\left(y_j , z_j\mid \theta\right)\right]=\log\left[P\left(y_1 , z_1\mid \theta\right)\right] +\sum_{j=2}^{n}\log\left[P\left(y_j , z_j\mid \theta\right)\right] j=1∑nlog[P(yj,zj∣θ)]=log[P(y1,z1∣θ)]+j=2∑nlog[P(yj,zj∣θ)]
则
∑
Z
∏
j
=
1
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
∑
j
=
1
n
log
[
P
(
y
j
,
z
j
∣
θ
)
]
=
∑
Z
{
log
[
P
(
y
1
,
z
1
∣
θ
)
]
+
∑
j
=
2
n
log
[
P
(
y
j
,
z
j
∣
θ
)
]
}
P
(
z
1
∣
y
1
,
θ
(
i
)
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
=
∑
Z
{
log
[
P
(
y
1
,
z
1
∣
θ
)
]
P
(
z
1
∣
y
1
,
θ
(
i
)
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
+
∑
j
=
2
n
log
[
P
(
y
j
,
z
j
∣
θ
)
]
P
(
z
1
∣
y
1
,
θ
(
i
)
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
}
\sum_{Z}\prod_{j=1}^{n} P\left(z_j|y_j,\theta^{(i)}\right) \sum_{j=1}^{n}\log\left[P\left(y_j , z_j\mid \theta\right)\right] =\\ \sum_{Z}\left\{\log\left[P\left(y_1 , z_1\mid \theta\right)\right] +\sum_{j=2}^{n}\log\left[P\left(y_j , z_j\mid \theta\right)\right]\right\}P\left(z_1|y_1,\theta^{(i)}\right)\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)\\ =\sum_{Z}\left\{\log\left[P\left(y_1 , z_1\mid \theta\right)\right] P\left(z_1|y_1,\theta^{(i)}\right)\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)+\\ \sum_{j=2}^{n}\log\left[P\left(y_j , z_j\mid \theta\right)\right] P\left(z_1|y_1,\theta^{(i)}\right)\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)\right\}
Z∑j=1∏nP(zj∣yj,θ(i))j=1∑nlog[P(yj,zj∣θ)]=Z∑{log[P(y1,z1∣θ)]+j=2∑nlog[P(yj,zj∣θ)]}P(z1∣y1,θ(i))j=2∏nP(zj∣yj,θ(i))=Z∑{log[P(y1,z1∣θ)]P(z1∣y1,θ(i))j=2∏nP(zj∣yj,θ(i))+j=2∑nlog[P(yj,zj∣θ)]P(z1∣y1,θ(i))j=2∏nP(zj∣yj,θ(i))}
哇槽,越弄公式越长,注意到
Z
Z
Z只有0和1两种取值,且括号内的两项为加的关系,可以直接分离,因此先看到第一项:
∑
Z
log
[
P
(
y
1
,
z
1
∣
θ
)
]
P
(
z
1
∣
y
1
,
θ
(
i
)
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
=
log
[
P
(
y
1
,
z
1
=
0
∣
θ
)
]
P
(
z
1
=
0
∣
y
1
,
θ
(
i
)
)
∑
Z
(
z
j
≠
1
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
+
log
[
P
(
y
1
,
z
1
=
1
∣
θ
)
]
P
(
z
1
=
1
∣
y
1
,
θ
(
i
)
)
∑
Z
(
z
j
≠
1
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
\sum_{Z}\log\left[P\left(y_1 , z_1\mid \theta\right)\right] P\left(z_1|y_1,\theta^{(i)}\right)\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)=\\ \log\left[P\left(y_1 , z_1=0\mid \theta\right)\right] P\left(z_1=0|y_1,\theta^{(i)}\right)\sum_{Z(z_j\neq1)}\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)+\\ \log\left[P\left(y_1 , z_1=1\mid \theta\right)\right] P\left(z_1=1|y_1,\theta^{(i)}\right)\sum_{Z(z_j\neq1)}\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)
Z∑log[P(y1,z1∣θ)]P(z1∣y1,θ(i))j=2∏nP(zj∣yj,θ(i))=log[P(y1,z1=0∣θ)]P(z1=0∣y1,θ(i))Z(zj=1)∑j=2∏nP(zj∣yj,θ(i))+log[P(y1,z1=1∣θ)]P(z1=1∣y1,θ(i))Z(zj=1)∑j=2∏nP(zj∣yj,θ(i))
这里出现了一个最振奋人心的项,
∑
Z
(
≠
1
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
\sum_{Z(\neq1)}\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)
∑Z(=1)∏j=2nP(zj∣yj,θ(i)),因为它等于1!!!于是有:
∑
Z
log
[
P
(
y
1
,
z
1
∣
θ
)
]
P
(
z
1
∣
y
1
,
θ
(
i
)
)
∏
j
=
2
n
P
(
z
j
∣
y
j
,
θ
(
i
)
)
=
log
[
P
(
y
1
,
z
1
=
0
∣
θ
)
]
P
(
z
1
=
0
∣
y
1
,
θ
(
i
)
)
+
log
[
P
(
y
1
,
z
1
=
1
∣
θ
)
]
P
(
z
1
=
1
∣
y
1
,
θ
(
i
)
)
=
∑
Z
log
[
P
(
y
1
,
z
1
∣
θ
)
]
P
(
z
1
∣
y
1
,
θ
(
i
)
)
\sum_{Z}\log\left[P\left(y_1 , z_1\mid \theta\right)\right] P\left(z_1|y_1,\theta^{(i)}\right)\prod_{j=2}^{n} P\left(z_j|y_j,\theta^{(i)}\right)=\\ \log\left[P\left(y_1 , z_1=0\mid \theta\right)\right] P\left(z_1=0|y_1,\theta^{(i)}\right)+\\ \log\left[P\left(y_1 , z_1=1\mid \theta\right)\right] P\left(z_1=1|y_1,\theta^{(i)}\right)\\= \sum_{Z}\log\left[P\left(y_1 , z_1\mid \theta\right)\right] P\left(z_1|y_1,\theta^{(i)}\right)
Z∑log[P(y1,z1∣θ)]P(z1∣y1,θ(i))j=2∏nP(zj∣yj,θ(i))=log[P(y1,z1=0∣θ)]P(z1=0∣y1,θ(i))+log[P(y1,z1=1∣θ)]P(z1=1∣y1,θ(i))=Z∑log[P(y1,z1∣θ)]P(z1∣y1,θ(i))
同理也可以求得第二项中各次的投币结果有着相同的形式,于是我们惊喜的得到:
Q
(
θ
,
θ
(
i
)
)
=
∑
j
=
1
n
∑
Z
log
[
P
(
y
j
,
z
j
∣
θ
)
]
P
(
z
j
∣
y
j
,
θ
(
i
)
)
Q\left(\theta, \theta^{(i)}\right)=\sum_{j=1}^{n}\sum_{Z}\log\left[P\left(y_j , z_j\mid \theta\right)\right] P\left(z_j|y_j,\theta^{(i)}\right)
Q(θ,θ(i))=j=1∑nZ∑log[P(yj,zj∣θ)]P(zj∣yj,θ(i))
将抛硬币模型代入其中有:
P
(
z
j
∣
y
j
,
θ
(
i
)
)
=
{
π
(
i
)
(
p
(
i
)
)
y
j
(
1
−
p
(
i
)
)
l
−
y
l
π
(
i
)
(
p
(
i
)
)
y
j
(
1
−
p
(
i
)
)
1
−
y
j
+
(
1
−
π
(
i
)
)
(
q
(
i
)
)
y
j
(
1
−
q
(
i
)
)
1
−
y
j
=
μ
j
(
i
+
1
)
,
z
j
=
1
1
−
μ
j
(
i
+
1
)
,
z
j
=
0
P\left(z_j|y_j,\theta^{(i)}\right) =\left\{\begin{array}{ll} \frac{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{l-y_{l}}}{\pi^{(i)}\left(p^{(i)}\right)^{y_{j}}\left(1-p^{(i)}\right)^{1-y_{j}}+\left(1-\pi^{(i)}\right)\left(q^{(i)}\right)^{y_{j}}\left(1-q^{(i)}\right)^{1-y_{j}}}=\mu_{j}^{(i+1)}&, z_{j}=1 \\ 1-\mu_{j}^{(i+1)},& z_{j}=0 \end{array}\right.
P(zj∣yj,θ(i))=⎩⎨⎧π(i)(p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))1−yjπ(i)(p(i))yj(1−p(i))l−yl=μj(i+1)1−μj(i+1),,zj=1zj=0
其中
μ
(
i
+
1
)
\mu^{(i+1)}
μ(i+1)为观测数据
y
j
y_j
yj来自抛硬币B的概率。
P
(
y
j
,
z
j
∣
θ
)
=
{
π
p
y
j
(
1
−
p
)
1
−
y
j
,
z
j
=
1
(
1
−
π
)
q
y
j
(
1
−
q
)
1
−
y
j
,
z
j
=
0
P\left(y_{j}, z_{j} \mid \theta\right)=\left\{\begin{array}{ll} \pi p^{y_{j}}(1-p)^{1-y_{j}}& , z_{j}=1 \\ (1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}, & z_{j}=0 \end{array}\right.
P(yj,zj∣θ)={πpyj(1−p)1−yj(1−π)qyj(1−q)1−yj,,zj=1zj=0
则
Q
(
θ
,
θ
(
i
)
)
=
∑
j
=
1
n
[
P
(
z
j
=
1
∣
y
j
,
θ
(
i
)
)
log
P
(
y
j
,
z
j
=
1
∣
θ
)
+
P
(
z
j
=
0
∣
y
j
,
θ
(
i
)
)
log
P
(
y
j
,
z
j
=
0
∣
θ
)
]
=
∑
j
=
1
n
{
μ
j
(
i
+
1
)
log
[
π
p
y
j
(
1
−
p
)
1
−
y
j
]
+
(
1
−
μ
j
(
i
+
1
)
)
log
[
(
1
−
π
)
q
y
j
(
1
−
q
)
1
−
y
j
]
}
Q\left(\theta, \theta^{(i)}\right)= \sum_{j=1}^{n}\left[P\left(z_j=1 \mid y_j, \theta^{(i)}\right) \log P(y_j, z_j=1 \mid \theta)+\\P\left(z_j=0 \mid y_j, \theta^{(i)}\right) \log P(y_j, z_j=0 \mid \theta)\right]\\ = \sum_{j=1}^{n}\left\{\mu_{j}^{(i+1)} \log [\pi p^{y_{j}}(1-p)^{1-y_{j}}]+ (1-\mu_{j}^{(i+1)}) \log[(1-\pi) q^{y_{j}}(1-q)^{1-y_{j}}]\right\}
Q(θ,θ(i))=j=1∑n[P(zj=1∣yj,θ(i))logP(yj,zj=1∣θ)+P(zj=0∣yj,θ(i))logP(yj,zj=0∣θ)]=j=1∑n{μj(i+1)log[πpyj(1−p)1−yj]+(1−μj(i+1))log[(1−π)qyj(1−q)1−yj]}
到此我们就完美的得到了一个相对比较简洁的
Q
Q
Q函数表达式,也就相当于得到了EM算法中E步的所有必要工具。
接下来的M步是求使得
Q
Q
Q函数极大化的
θ
\theta
θ,事实上就是求各参数对
Q
Q
Q函数的偏导数为0的点:
∂
Q
∂
π
=
∑
j
=
1
n
[
μ
j
(
i
+
1
)
π
−
1
−
μ
j
(
i
+
1
)
1
−
π
]
=
∑
j
=
1
n
[
μ
j
(
i
+
1
)
−
μ
j
(
i
+
1
)
π
π
(
1
−
π
)
−
π
−
μ
j
(
i
+
1
)
π
π
(
1
−
π
)
]
=
∑
j
=
1
n
[
μ
j
(
i
+
1
)
−
π
π
(
1
−
π
)
]
=
∑
j
=
1
n
μ
j
(
i
+
1
)
−
n
π
π
(
1
−
π
)
=
0
\frac{\partial Q}{\partial \pi}=\sum_{j=1}^{n}\left[\frac{\mu_{j}^{(i+1)} }{\pi}- \frac{1-\mu_{j}^{(i+1)} }{ 1-\pi }\right]= \sum_{j=1}^{n}\left[\frac{\mu_{j}^{(i+1)}-\mu_{j}^{(i+1)}\pi }{\pi(1-\pi)}- \frac{\pi-\mu_{j}^{(i+1)}\pi }{ \pi(1-\pi) }\right]\\= \sum_{j=1}^{n}\left[\frac{\mu_{j}^{(i+1)}-\pi }{\pi(1-\pi)}\right]= \frac{\sum_{j=1}^{n}\mu_{j}^{(i+1)}-n\pi }{\pi(1-\pi)}=0
∂π∂Q=j=1∑n[πμj(i+1)−1−π1−μj(i+1)]=j=1∑n[π(1−π)μj(i+1)−μj(i+1)π−π(1−π)π−μj(i+1)π]=j=1∑n[π(1−π)μj(i+1)−π]=π(1−π)∑j=1nμj(i+1)−nπ=0
有
π
i
+
1
=
1
n
∑
j
=
1
n
μ
j
(
i
+
1
)
\pi^{i+1}=\frac{1}{n}\sum_{j=1}^{n}\mu_{j}^{(i+1)}
πi+1=n1j=1∑nμj(i+1)
同理可得:
p
(
i
+
1
)
=
∑
j
=
1
n
μ
j
(
i
+
1
)
y
j
∑
j
=
1
n
μ
j
(
i
+
1
)
p^{(i+1)}=\frac{\sum_{j=1}^{n} \mu_{j}^{(i+1)} y_{j}}{\sum_{j=1}^{n} \mu_{j}^{(i+1)}}
p(i+1)=∑j=1nμj(i+1)∑j=1nμj(i+1)yj
q ( i + 1 ) = ∑ j = 1 n ( 1 − μ j ( i + 1 ) ) y j ∑ j = 1 n ( 1 − μ j ( i + 1 ) ) q^{(i+1)}=\frac{\sum_{j=1}^{n}\left(1-\mu_{j}^{(i+1)}\right) y_{j}}{\sum_{j=1}^{n}\left(1-\mu_{j}^{(i+1)}\right)} q(i+1)=∑j=1n(1−μj(i+1))∑j=1n(1−μj(i+1))yj
到此,M步所需要的所有东西也都齐全了,问题解决!
EM算法解决抛硬币问题—实现代码
如果是熟悉python编程的老哥可以去下载黄海广博士关于《统计学习方法》这本书的实现代码,非常给力。下面代码是本人使用C++写的:
#include <iostream>
#include <cmath>
#define N 10
//设置变量并初始化
//double pi = 0.5, p = 0.5, q = 0.5;
double pi = 0.4, p = 0.6, q = 0.7;
int y[N] = { 1,1,0,1,0,0,1,0,1,1 };
double u[N];
void EStep()
{
for (size_t i = 0; i < N; i++)
{
u[i] = pi * std::pow(p, y[i]) * std::pow(1 - p, 1 - y[i]);
u[i] = u[i] / (u[i] + (1 - pi) * std::pow(q, y[i]) * std::pow(1 - q, 1 - y[i]));
}
}
void MStep()
{
double addu = 0;
for (size_t i = 0; i < N; i++)
{
addu += u[i];
}
pi = 1 / (double)N * addu;
p = 0;
for (size_t i = 0; i < N; i++)
{
p+= u[i]*y[i];
}
p = p / addu;
double adduMOne = 0;
for (size_t i = 0; i < N; i++)
{
adduMOne += 1-u[i];
}
q = 0;
for (size_t i = 0; i < N; i++)
{
q += (1-u[i]) * y[i];
}
q = q / adduMOne;
}
int main()
{
std::cout << "原始数据" << std::endl;
std::cout << "pi=" << pi << " " << "p=" << p << " " << "q=" << q << std::endl << std::endl;
EStep();
MStep();
std::cout << "第一次迭代结果" << std::endl;
std::cout << "pi=" << pi << " " << "p=" << p << " " << "q=" << q << std::endl << std::endl;
EStep();
MStep();
std::cout << "第二次迭代结果" << std::endl;
std::cout << "pi=" << pi << " " << "p=" << p << " " << "q=" << q << std::endl << std::endl;
}
扩展知识—符号/名词含义
P ( y ∣ θ ) P(y \mid \theta) P(y∣θ)->在 θ \theta θ发生的情况下, y y y发生的概率
P ( y , z ∣ θ ) P(y, z\mid\theta) P(y,z∣θ)->在 θ \theta θ发生的情况下, y y y和 z z z发生的概率
P ( y ∣ z , θ ) P(y \mid z, \theta) P(y∣z,θ)->在 z z z和 θ \theta θ发生的情况下, y y y发生的概率
解析解->指通过严格的公式所求得的解,当给出解的具体函数形式,从解的表达式中就可以算出任何对应值。
比如等式:x^2=2,其解析解为:x=sqrt(2),而数值解为:x=1.414
参考文章
https://zhuanlan.zhihu.com/p/133752509
https://blog.csdn.net/zsdust/article/details/100042491?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param