EM算法是一种迭代算法,用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。
EM算法每次迭代由两步组成:E步,求期望,M步,求极大。称为期望极大算法(expectation maximization)
观测数据为
Y
=
(
Y
,
Y
2
,
…
,
Y
n
)
T
Y=(Y_,Y_2,\dots,Y_n)^T
Y=(Y,Y2,…,Yn)T,不可观测数据为
Z
=
(
Z
1
,
Z
2
,
…
,
Z
n
)
T
Z=(Z_1,Z_2,\dots,Z_n)^T
Z=(Z1,Z2,…,Zn)T,则观测数据的似然函数为
P
(
Y
∣
θ
)
=
∑
z
P
(
Z
∣
θ
)
P
(
Y
∣
Z
,
θ
)
P(Y|\theta)=\sum_zP(Z|\theta)P(Y|Z,\theta)
P(Y∣θ)=z∑P(Z∣θ)P(Y∣Z,θ)
=
∏
j
=
n
[
π
p
y
j
(
1
−
p
)
1
−
y
i
+
(
1
−
π
)
q
y
j
(
1
−
q
)
1
−
y
j
]
=\prod_{j=}^n[\pi p^{y_j}(1-p)^{1-y_i}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]
=j=∏n[πpyj(1−p)1−yi+(1−π)qyj(1−q)1−yj]
模型参数
θ
=
(
π
,
p
,
q
)
\theta=(\pi,p,q)
θ=(π,p,q)的极大似然估计为
θ
^
=
a
r
g
m
a
x
θ
l
o
g
P
(
Y
∣
θ
)
\hat\theta=argmax_\theta logP(Y|\theta)
θ^=argmaxθlogP(Y∣θ)
这个问题没有解析解,只有通过迭代的方法求解。
EM算法首先选取参数的初值,记做
θ
(
0
)
=
(
π
(
0
)
,
p
(
0
)
,
q
(
0
)
)
\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})
θ(0)=(π(0),p(0),q(0)),然后通过下面的步骤接待计算参数的估计值,直至收敛为止。第i次迭代参数的估计值为
θ
(
i
)
=
(
π
(
i
)
,
p
(
i
)
,
q
(
i
)
)
\theta^{(i)}=(\pi^{(i)},p^{(i)},q^{(i)})
θ(i)=(π(i),p(i),q(i))。EM算法第i+1次迭代如下:
E步(求期望):计算在模型参数
π
(
i
)
,
p
(
i
)
,
q
(
i
)
\pi^{(i)},p^{(i)},q^{(i)}
π(i),p(i),q(i)下观测数据
y
j
y_j
yj来自掷硬币B的概率
μ
j
(
i
+
1
)
=
π
(
i
)
(
p
(
i
)
)
y
j
(
1
−
p
(
i
)
)
1
−
y
j
π
(
i
)
(
p
(
i
)
)
y
j
(
1
−
p
(
i
)
)
1
−
y
j
+
(
1
−
π
(
i
)
)
(
q
(
i
)
)
y
j
(
1
−
q
(
i
)
)
1
−
y
j
\mu_j^{(i+1)}=\frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}+(1-\pi^{(i)})(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}}
μj(i+1)=π(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))1−yj
M步(求极大化):计算模型参数的新估计值
π
(
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μji+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(1-\mu_j^{(i+1)})y_j}{\sum_{j=1}^n(1-\mu_j^{(i+1)})}
q(i+1)=∑j=1n(1−μj(i+1))∑j=1n(1−μj(i+1))yj
EM算法
输入:观测变量数据Y,隐变量数据Z,联合分布
P
(
Y
,
Z
∣
θ
)
P(Y,Z|\theta)
P(Y,Z∣θ),条件分布
P
(
Z
∣
Y
,
θ
)
P(Z|Y,\theta)
P(Z∣Y,θ)
输出:模型参数θ
(1)选择参数的初值
θ
(
0
)
\theta^{(0)}
θ(0),开始迭代
(2)E步:记
θ
(
i
)
\theta^{(i)}
θ(i)为第i 次迭代参数θ的估计值,在第i+1次迭代的E步,计算
Q
(
θ
,
θ
(
i
)
)
=
E
Z
[
l
o
g
P
(
Y
,
Z
∣
θ
)
∣
Y
,
θ
(
i
)
]
=
∑
z
l
o
g
P
(
Y
,
Z
∣
θ
)
P
(
Z
∣
Y
,
θ
(
i
)
)
Q(\theta,\theta^{(i)})=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}]=\sum_zlogP(Y,Z|\theta)P(Z|Y,\theta^{(i)})
Q(θ,θ(i))=EZ[logP(Y,Z∣θ)∣Y,θ(i)]=z∑logP(Y,Z∣θ)P(Z∣Y,θ(i))
这里
P
(
Z
∣
Y
,
θ
(
i
)
)
P(Z|Y,\theta^{(i)})
P(Z∣Y,θ(i))是在给定观测数据Y和当前参数估计
θ
(
i
)
\theta^{(i)}
θ(i)下隐变量数据Z的条件概率分布;
(3)M步:求使
Q
(
θ
,
θ
(
i
)
)
Q(\theta,\theta^{(i)})
Q(θ,θ(i))极大化的
θ
\theta
θ,确定第i+1次迭代的参数估计值
θ
(
i
+
1
)
\theta^{(i+1)}
θ(i+1),
θ
(
i
+
1
)
=
a
r
g
m
a
x
θ
Q
(
θ
,
θ
(
i
)
)
\theta^{(i+1)}=argmax_\theta Q(\theta,\theta^{(i)})
θ(i+1)=argmaxθQ(θ,θ(i))
(4)重复第(2)(3)步直至收敛。迭代停止条件,一般是对较小的正整数
ε
1
,
ε
2
\varepsilon_1,\varepsilon_2
ε1,ε2满足
∣
∣
θ
(
i
+
1
)
−
θ
(
i
)
<
ε
1
||\theta^{(i+1)}-\theta^{(i)}<\varepsilon_1
∣∣θ(i+1)−θ(i)<ε1或
∣
∣
Q
(
θ
(
i
+
1
)
,
θ
(
i
)
)
−
Q
(
θ
(
i
)
,
θ
(
i
)
)
∣
∣
<
ε
2
||Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})||<\varepsilon_2
∣∣Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∣∣<ε2,则停止迭代。
EM算法在高斯混合模型学习中的应用
高斯混合模型
定义:高斯混合模型具有如下形式的概率分布:
P
(
y
∣
θ
)
=
∑
k
=
1
K
α
k
ψ
(
y
∣
θ
k
)
P(y|\theta)=\sum_{k=1}^K\alpha_k\psi(y|\theta_k)
P(y∣θ)=∑k=1Kαkψ(y∣θk),其中
α
k
\alpha_k
αk是系数,
α
k
≥
0
,
∑
k
=
1
K
α
k
=
1
;
ψ
(
y
∣
θ
k
)
\alpha_k\geq0,\sum_{k=1}^K\alpha_k=1;\psi(y|\theta_k)
αk≥0,∑k=1Kαk=1;ψ(y∣θk)是高斯分布密度,
θ
k
=
(
μ
k
,
σ
k
2
)
\theta_k=(\mu_k,\sigma_k^2)
θk=(μk,σk2),
ψ
(
y
∣
θ
k
)
=
1
2
π
σ
k
e
x
p
(
−
(
y
−
μ
k
)
2
2
σ
k
2
)
\psi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}exp\bigg(-\frac{(y-\mu_k)^2}{2\sigma_k^2}\bigg)
ψ(y∣θk)=2πσk1exp(−2σk2(y−μk)2)称为第k个分模型。
算法:高斯混合模型参数估计的EM算法
输入:观测数据
y
1
,
y
2
,
…
,
y
N
y_1,y_2,\dots,y_N
y1,y2,…,yN,高斯混合模型;
输出:高斯混合模型参数
(1)取参数的初始值开始迭代
(2)E步:依据当前模型参数,计算分模型k对观测数据
y
j
y_j
yj的响应度
γ
^
j
k
=
α
k
ψ
(
y
j
∣
θ
k
)
∑
k
=
1
K
α
k
ψ
(
y
j
∣
θ
k
)
,
j
=
1
,
2
,
…
,
N
;
k
=
1
,
2
,
…
,
K
\hat\gamma_{jk}=\frac{\alpha_k\psi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\psi(y_j|\theta_k)},j=1,2,\dots,N;k=1,2,\dots,K
γ^jk=∑k=1Kαkψ(yj∣θk)αkψ(yj∣θk),j=1,2,…,N;k=1,2,…,K
(3)M步:计算新一轮迭代的模型参数
μ
^
k
=
∑
j
=
1
N
γ
^
j
k
y
j
∑
j
=
1
N
γ
^
j
k
,
k
=
1
,
2
,
…
,
K
\hat\mu_k=\frac{\sum_{j=1}^N\hat\gamma_{jk}y_j}{\sum_{j=1}^N\hat\gamma_{jk}},k=1,2,\dots,K
μ^k=∑j=1Nγ^jk∑j=1Nγ^jkyj,k=1,2,…,K
σ
^
k
2
=
∑
j
=
1
N
γ
^
j
k
(
y
j
−
μ
k
)
2
∑
j
=
1
N
γ
^
j
k
,
k
=
1
,
2
,
…
,
K
\hat\sigma_k^2=\frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}},k=1,2,\dots,K
σ^k2=∑j=1Nγ^jk∑j=1Nγ^jk(yj−μk)2,k=1,2,…,K
α
^
k
=
∑
j
=
1
N
γ
^
j
k
N
,
k
=
1
,
2
,
…
,
K
\hat\alpha_k=\frac{\sum_{j=1}^N\hat\gamma_{jk}}{N},k=1,2,\dots,K
α^k=N∑j=1Nγ^jk,k=1,2,…,K
(4)重复第(2)(3)步直至收敛。