来自https://www.jianshu.com/p/f4ff48025c52
https://blog.csdn.net/wendaomudong_l2d4/article/details/79005461
1、引言
概率模型有时既含有观测变量,又含有隐变量或潜在变量,如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计法估计模型参数,但是,当模型中含有隐变量时,就不能简单的使用这些方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法。
EM算法流程
E步:对完全数据的对数似然函数$\log(P(Y,Z|\theta)$
求关于$P(Z|Y,\theta^{(i)})$
的数学期望,
E
Z
∣
Y
,
θ
(
i
)
[
l
o
g
(
P
(
Y
,
Z
∣
θ
)
)
]
\quad E_{Z|Y,\theta^{(i)}}[log(P(Y,Z|\theta))]
EZ∣Y,θ(i)[log(P(Y,Z∣θ))]
其中$\theta^{(i)}$
是第i次迭代时,$
θ
\theta
θ`的估计值
M步:对E步的结果求极值
EM算法步骤
一般的,用Y表示观测随机变量的数据,Z表示隐随机变量的数据,Y和Z连在一起称为完全数据,只有观测数据Y称为不完全数据,假设给定观测数据Y,其概率分布为P(Y|θ),那么不完全数据的似然函数就是P(Y|θ),对数似然函数是L(θ) = log(P(Y|θ)),假设Y和Z的联合概率分布是P(Y,Z|θ),那么完全数据的对数似然函数是logP(Y,Z|θ)。
EM算法通过迭代求L(θ) = log(P(Y|θ))的极大似然估计,每次迭代包括两步:E步,求期望,M步,求最大化,下面介绍EM算法的步骤:
KaTeX parse error: Can't use function '$' in math mode at position 6: \quad$̲`输入:观测变量数据Y,隐变量…,条件分布` P ( Z ∣ Y , θ ) P(Z|Y,\theta) P(Z∣Y,θ)$;
$
\quad
输出:模型参数
θ
\theta
θ`.
KaTeX parse error: Can't use function '$' in math mode at position 10: \quad (1)$̲`选择参数的初值`$\thet…,开始迭代;
$
(
2
)
\quad (2)
(2)E步:记
θ
(
i
)
\theta^{(i)}
θ(i)为第i次迭代参数
θ
\theta
θ`的估计值,在第i+1次迭代的E步,计算
Q ( θ , θ ( i ) ) = E Y , θ ( i ) [ log ( P ( Y , Z ∣ θ ) ) ] \qquad\qquad Q(\theta,\theta^{(i)})=E_{Y,\theta^{(i)}}[\log( P(Y,Z|\theta))] Q(θ,θ(i))=EY,θ(i)[log(P(Y,Z∣θ))]
= ∑ Z l o g ( P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) ) \qquad\qquad\qquad\quad\ =\sum\limits_Z log(P(Y,Z|\theta)P(Z|Y,\theta^{(i)})) =Z∑log(P(Y,Z∣θ)P(Z∣Y,θ(i)))
这里,$P(Z|Y,\theta^{(i)})$
是在给定观测数据Y和当前的参数估计$\theta^{(i)}$
下隐变量数据Z的条件概率分布;
KaTeX parse error: Can't use function '$' in math mode at position 6: \quad$̲`(3) M步: 求使`$Q(…
θ ( i + 1 ) = arg max θ Q ( θ , θ ( i ) ) \qquad\qquad\qquad \theta^{(i+1)}=\arg\underset{\theta}{\max}Q(\theta,\theta^{(i)}) θ(i+1)=argθmaxQ(θ,θ(i))
$ \quad `(4) 重复第(2)步和第(3)步,直到收敛.
2、三硬币模型描述
三硬币问题是这样的:
假设有三枚硬币,分别记为A、B、C。这些硬币正面的概率分别为π,p,q,进行如下的抛硬币实验:先掷硬币A,根据其结果选出硬币B或者硬币C,正面选硬币B,反面选硬币C,然后掷选出的硬币,掷硬币的记录,出现正面记作1,出现反面记作0,独立地重复n次实验(这里n=10),然后观测结果如下:
1,1,0,1,0,0,1,0,1,1
假设只能观测到掷硬币的结果,不能观测掷硬币的过程,问如何估计三硬币正面出现的概率,即三硬币模型的参数π,p,q。
3、三硬币问题表示
符号标记:
KaTeX parse error: Can't use function '$' in math mode at position 10: \quad y_j$̲`为第j次实验的观测 \quad Z$为隐变量,表示掷硬币A出现的结果.该变量只有两个取值0,1 $$\quad z_j$
为第j次实验时,掷硬币A出现的结果,$z_j=1$
表示硬币A掷出正面,$z_j=0$
表示硬币A掷出反面
KaTeX parse error: Can't use function '$' in math mode at position 12: \quad\theta$̲`表示参数集合`$\pi,p,…\quad\theta^{(i)}
‘
为
第
i
次
迭
代
时
,
`为第i次迭代时,
‘为第i次迭代时,
π
,
p
,
q
\pi,p,q
π,p,q`的估计值
4、EM算法求解三硬币模型
假设当前模型的参数为π,p,q时,隐含变量来自于硬币A的后验概率,
对数似然函数为:
log ( P ( Y , Z ∣ θ ) ) = log ( ∏ j = 1 n p ( y j , z j ∣ θ ) ) = ∑ j = 1 n log ( p ( y j , z j ∣ θ ) ) \log(P(Y,Z|\theta))=\log(\prod\limits_{j=1}^n p(y_j,z_j|\theta))=\sum\limits_{j=1}^n \log(p(y_j,z_j|\theta)) log(P(Y,Z∣θ))=log(j=1∏np(yj,zj∣θ))=j=1∑nlog(p(yj,zj∣θ))
期望为:
` E Z ∣ Y , θ ( i ) [ log ( P ( Y , Z ∣ θ ) ) ] E_{Z|Y,\theta^{(i)}}[\log(P(Y,Z|\theta))] EZ∣Y,θ(i)[log(P(Y,Z∣θ))]$
= ∑ j = 1 n ∑ z j [ p ( z j ∣ y j , θ ( i ) ) log ( p ( y j , z j ∣ θ ) ) ] =\sum\limits_{j=1}^n\sum\limits_{z_j}[p(z_j|y_j,\theta^{(i)})\log(p(y_j,z_j|\theta))] =j=1∑nzj∑[p(zj∣yj,θ(i))log(p(yj,zj∣θ))]
= ∑ 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 ∣ θ ) ) ) =\sum\limits_{j=1}^n(p(z_j=1|y_j,\theta^{(i)})\log(p(y_j,z_j=1|\theta))+p(z_j=0|y_j,\theta^{(i)})\log(p(y_j,z_j=0|\theta))) =j=1∑n(p(zj=1∣yj,θ(i))log(p(yj,zj=1∣θ))+p(zj=0∣yj,θ(i))log(p(yj,zj=0∣θ)))
对于概率`
p
(
z
j
,
y
j
∣
θ
)
p(z_j,y_j|\theta)
p(zj,yj∣θ)$:
p
(
y
j
,
z
j
=
1
∣
θ
)
=
p
(
y
j
∣
z
j
=
1
,
θ
)
p
(
z
j
=
1
∣
θ
)
=
π
p
y
j
(
1
−
p
)
1
−
y
j
\quad p(y_j,z_j=1|\theta)=p(y_j|z_j=1,\theta)p(z_j=1|\theta)=\pi p^{y_j}(1-p)^{1-y_j}
p(yj,zj=1∣θ)=p(yj∣zj=1,θ)p(zj=1∣θ)=πpyj(1−p)1−yj,
p ( y j , z j = 0 ∣ θ ) = p ( y j ∣ z j = 0 , θ ) p ( z j = 0 ∣ θ ) = ( 1 − π ) q y j ( 1 − q ) 1 − y j \quad p(y_j,z_j=0|\theta)=p(y_j|z_j=0,\theta)p(z_j=0|\theta)=(1-\pi) q^{y_j}(1-q)^{1-y_j} p(yj,zj=0∣θ)=p(yj∣zj=0,θ)p(zj=0∣θ)=(1−π)qyj(1−q)1−yj
所以最终结果:
Q ( θ , θ ( i ) ) = E Z ∣ Y , θ ( i ) [ log ( P ( Y , Z ∣ θ ) ) ] \qquad Q(\theta,\theta^{(i)})=E_{Z|Y,\theta^{(i)}}[\log(P(Y,Z|\theta))] Q(θ,θ(i))=EZ∣Y,θ(i)[log(P(Y,Z∣θ))]
= ∑ j = 1 n ( μ j ( i ) ∗ log ( π p y j ( 1 − p ) 1 − y j ) + ( 1 − μ j ( i ) ) ⋆ log ( ( 1 − π ) q y j ( 1 − q ) 1 − y j ) ) \qquad\qquad = \sum\limits_{j=1}^n(\mu_j^{(i)} * \log(\pi p^{y_j}(1-p)^{1-y_j})+(1-\mu_j^{(i)})\star \log((1-\pi)q^{y_j}(1-q)^{1-y_j})) =j=1∑n(μj(i)∗log(πpyj(1−p)1−yj)+(1−μj(i))⋆log((1−π)qyj(1−q)1−yj))
上式中` μ j ( i ) = p ( z j = 1 ∣ y j ; θ ( i ) ) \mu_j^{(i)}=p(z_j=1|y_j;\theta^{(i)}) μj(i)=p(zj=1∣yj;θ(i))$.
EM算法的具体实现
E-Step:
μ j ( i ) = π p y j ( 1 − p ) 1 − y j π p y j ( 1 − p ) 1 − y j + ( 1 − π ) q y j ( 1 − q ) 1 − y j \qquad \mu_j^{(i)}=\frac{\pi p^{y_j}(1-p)^{1-y_j}}{\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}} μj(i)=πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yjπpyj(1−p)1−yj,
那么隐含变量来自于硬币C的后验概率自然为` 1 − μ j ( i ) 1-\mu_j^{(i)} 1−μj(i)$.
def cal_u(phi, p, q, yj):
return phi * math.pow(p, yj) * math.pow(1 - p, 1 - yj) / \
float(phi * math.pow(p, yj) * math.pow(1 - p, 1 - yj) +
(1 - phi) * math.pow(q, yj) * math.pow(1 - q, 1 - yj))
def e_step(phi,p,q,y):
"""
e步计算
:param phi: 下一次迭代开始的 pi
:param p: 下一次迭代开始的 p
:param q: 下一次迭代开始的 q
:param y: 观察数据
:return:
"""
return [cal_u(phi,p,q,yj) for yj in y]
M-Step:
下面分别对$\pi,p,q$
求偏导:
(1) 估计参数$\pi$
$$\qquad ‘ 对 ‘ `对` ‘对‘\pi$`求偏导:
\begin{array}{lll}
\frac{\partial Q}{\partial \pi}&=&\sum\limits_{j=1}^n(\mu_j^{(i)}* \frac{1}{\pi}-(1-\mu_j^{(i)})* \frac{1}{1-\pi})\\
&=&\sum\limits_{j=1}^n(\frac{\pi-\mu_j^{(i)}}{\pi(1-\pi)})\\
&=&\frac{n\pi-\sum\limits_{j=1}^n\mu_j^{(i)}}{\pi(1-\pi)}=0
\end{array}
所以π的估计为:
π = 1 n ∑ j = 1 n μ j ( i ) \qquad\qquad \pi=\frac{1}{n}\sum\limits_{j=1}^n \mu_j^{(i)} π=n1j=1∑nμj(i)
(2) 估计参数p
对p求偏导:
\begin{array}{lll}
\frac{\partial Q}{\partial p}&=&\sum\limits_{j=1}^n\mu_j^{(i)}*\frac{\pi(y_j p^{y_j-1}(1-p)^{1-y_j}+p^{y_j}(-(1-y_j)(1-p)^{-y_j})}{\pi p^{y_j}(1-p)^{1-y_j}}\\
&=&\sum\limits_{j=1}^n \mu_j^{(i)}* \frac{y_j p^{y_j-1}(1-p)^{-y_j}(1-p)+p^{y_j-1}*p(y_j-1)(1-p)^{-y_j}}{p^{y_j}(1-p)^{1-y_j}}\\
&=&\sum\limits_{j=1}^n\mu_j^{(i)}*\frac{y_j(1-p)+p*(y_j-1)}{p(1-p)}\\
&=&\sum\limits_{j=1}^n\mu_j^{(i)}*\frac{y_j-p}{p(1-p)}=0
\end{array}
p的估计为:
p = ∑ j = 1 n μ j ( i ) y j ∑ j = 1 n μ j ( i ) \qquad p=\frac{\sum\limits_{j=1}^n\mu_j^{(i)}y_j}{\sum\limits_{j=1}^n\mu_j^{(i)}} p=j=1∑nμj(i)j=1∑nμj(i)yj
(3) 估计参数q
对q求偏导:
\begin{array}{lll}
\frac{\partial Q}{\partial q}&=&\sum\limits_{j=1}^n(1-\mu_j^{(i)})*\frac{(1-\pi)(y_j q^{y_j-1}(1-q)^{1-y_j}+q_{y_j}(-(1-y_j)(1-q){-y_j}))}{(1-\pi)q^{y_j}(1-q)^{1-y_j}}\\
&=&\sum\limits_{j=1}^n(1-\mu_j^{(i)})*\frac{y_j q^{y_j-1}(1-q)^{-y_j}*(1-q)+q^{y_j-1}*q(y_j-1)(1-q)^{-y_j}}{q^{y_j}(1-q)^{1-y_j}}\\
&=&\sum\limits_{j=1}^n(1-\mu_j^{(i)})*\frac{y_j(1-q)+q*(y_j-1)}{q(1-q)}\\
&=&\sum\limits_{j=1}^n(1-\mu_j^{(i)})*\frac{y_j-q}{p(1-q)}=0
\end{array}
q的估计为:
p = ∑ j = 1 n ( 1 − μ j ( i ) ) y j ∑ j = 1 n ( 1 − μ j ( i ) ) \qquad\ p=\frac{\sum\limits_{j=1}^n(1-\mu_j^{(i)})y_j}{\sum\limits_{j=1}^n(1-\mu_j^{(i)})} p=j=1∑n(1−μj(i))j=1∑n(1−μj(i))yj
更新$\pi,p,q$
的代码如下:
def m_step(u,y):
phi1=sum(u)/len(u)
p1=sum([u[i]*y[i] for i in range(len(u))]) / sum(u)
q1=sum([(1-u[i])*y[i] for i in range(len(u))]) / sum([1-u[i] for i in range(len(u))])
return [phi1,p1,q1]
算法首先选取参数初始值θ(0)=π(0),p(0),q(0) ,然后迭代到收敛为止.
完整的python代码如下:
# -*- coding: utf-8 -*-
import math
def cal_u(phi, p, q, yj):
"""
u值计算
:param phi: 下一次迭代开始的 pi
:param p: 下一次迭代开始的 p
:param q: 下一次迭代开始的 q
:param yj: 观察数据第j个值,从0开始
:return:
"""
return phi * math.pow(p, yj) * math.pow(1 - p, 1 - yj) / \
float(phi * math.pow(p, yj) * math.pow(1 - p, 1 - yj) +
(1 - phi) * math.pow(q, yj) * math.pow(1 - q, 1 - yj))
def e_step(phi,p,q,y):
"""
e步计算
:param phi: 下一次迭代开始的 pi
:param p: 下一次迭代开始的 p
:param q: 下一次迭代开始的 q
:param y: 观察数据
:return:
"""
return [cal_u(phi,p,q,yj) for yj in y]
def m_step(u,y):
"""
m步计算
:param u: m步计算的u
:param y: 观察数据
:return:
"""
phi1=sum(u)/len(u)
p1=sum([u[i]*y[i] for i in range(len(u))]) / sum(u)
q1=sum([(1-u[i])*y[i] for i in range(len(u))]) / sum([1-u[i] for i in range(len(u))])
return [phi1,p1,q1]
def run(observed_y, start_phi, start_p, start_q, iter_num):
"""
:param observed_y: 观察数据
:param start_phi: 下一次迭代开始的pi $\pi$
:param start_p: 下一次迭代开始的p
:param start_q: 下一次迭代开始的q
:param iter_num: 迭代次数
:return:
"""
for i in range(iter_num):
u=e_step(start_phi, start_p, start_q, observed_y)
print (i,[start_phi,start_p,start_q])
if [start_phi,start_p,start_q]==m_step(u, observed_y):
break
else:
[start_phi,start_p,start_q]=m_step(u, observed_y)
if __name__ =="__main__":
# 观察数据
y = [1, 1, 0, 1, 0, 0, 1, 0, 1, 1]
# 初始化 pi,p q
[phi, p, q] = [0.4, 0.6, 0.7]
# 迭代计算
run(y,phi,p,q,100)
结果如下:
0 [0.4, 0.6, 0.7]
1 [0.40641711229946526, 0.5368421052631579, 0.6432432432432431]
2 [0.40641711229946537, 0.5368421052631579, 0.6432432432432431]
3 [0.40641711229946537, 0.536842105263158, 0.6432432432432431]