概述
文章目录
- 参考资料
- 1. 最大似然估计
- 1.1 原理
- 1.2 示例
- 2. EM算法
- 2.1 原理
- 2.2 示例
参考资料
- 统计计算中的优化问题
1. 最大似然估计
1.1 原理
统计中许多问题的计算最终都归结为一个最优化问题, 典型代表是最大似然估计(MLE)、各种拟似然估计方法、 非线性回归、惩罚函数方法(如svm、lasso)等。
最大似然估计经常需要用最优化算法计算, 最大似然估计问题有自身的特点, 可以直接用一般优化方法进行最大似然估计的计算, 但是利用最大似然估计的特点可以得到更有效的算法。
设总体
X
boldsymbol{X}
X 有概率密度(连续型随机变量)或概率分布(离散型随机变量)
p
(
x
∣
θ
)
,
θ
p(boldsymbol{x} mid boldsymbol{theta}), boldsymbol{theta}
p(x∣θ),θ 为
m
m
m 维的分布参数。有了一组样本
X
1
,
X
2
,
…
,
X
n
boldsymbol{X}_1, boldsymbol{X}_2, ldots, boldsymbol{X}_n
X1,X2,…,Xn 后,似然函数 为
L
(
θ
)
=
∏
i
=
1
n
p
(
X
i
∣
θ
)
L(boldsymbol{theta})=prod_{i=1}^n pleft(boldsymbol{X}_i mid boldsymbol{theta}right)
L(θ)=i=1∏np(Xi∣θ)
对数似然函数为
l
(
θ
)
=
∑
i
=
1
n
ln
p
(
X
i
∣
θ
)
l(boldsymbol{theta})=sum_{i=1}^n ln pleft(boldsymbol{X}_i mid boldsymbol{theta}right)
l(θ)=i=1∑nlnp(Xi∣θ)
最大化的步骤通过对 l ( θ ) l(boldsymbol{theta}) l(θ)求导等于0来解得。
1.2 示例
-
例1:设 X ∼ b ( 1 , p ) , X 1 , … , X n X sim b(1, p), X_1, ldots, X_n X∼b(1,p),X1,…,Xn 是来自 X X X 的一个样本, x 1 , … , x n x_1, ldots, x_n x1,…,xn 为观察值。试求参数 p p p 的最大似然 估计。
解:可知 X X X 的分布律为:
P { X = x } = p x ( 1 − p ) 1 − x P{X=x}=p^x(1-p)^{1-x} P{X=x}=px(1−p)1−x
故似然函数为:
L ( p ∣ x ) = ∏ i = 1 n p x i ( 1 − p ) 1 − x i = p ∑ i = 1 n x i ( 1 − p ) n − ∑ i = 1 n x i ln L ( p ∣ x ) = ∑ i = 1 n x i ln p + ( n − ∑ i = 1 n x i ) ln ( 1 − p ) begin{aligned} &L(p mid boldsymbol{x})=prod_{i=1}^n p^{x_i}(1-p)^{1-x_i}=p^{sum_{i=1}^n x_i}(1-p)^{n-sum_{i=1}^n x_i} \ &ln L(p mid boldsymbol{x})=sum_{i=1}^n x_i ln p+left(n-sum_{i=1}^n x_iright) ln (1-mathrm{p}) end{aligned} L(p∣x)=i=1∏npxi(1−p)1−xi=p∑i=1nxi(1−p)n−∑i=1nxilnL(p∣x)=i=1∑nxilnp+(n−i=1∑nxi)ln(1−p)
对待估参数求导为 0 , 有:
d d p ln L ( p ∣ x ) = ∑ i = 1 n x i p − n − ∑ i = 1 n x i 1 − p = 0 frac{d}{d p} ln L(p mid boldsymbol{x})=frac{sum_{i=1}^n x_i}{p}-frac{n-sum_{i=1}^n x_i}{1-p}=0 dpdlnL(p∣x)=p∑i=1nxi−1−pn−∑i=1nxi=0
解得 p p p 的最大似然估计值为 p ^ = x ˉ hat{p}=bar{x} p^=xˉ ,最大估计量为 p ^ = X ˉ hat{p}=bar{X} p^=Xˉ. -
例2: 设 X ∼ N ( μ , σ 2 ) , μ , σ 2 X sim Nleft(mu, sigma^2right), mu, sigma^2 X∼N(μ,σ2),μ,σ2 为未知参数, x 1 , … , x n x_1, ldots, x_n x1,…,xn 为来自 X X X 的一组观察值, 求 μ , σ 2 mu, sigma^2 μ,σ2 的最大似然估计量。
解: X X X 的概率密度为
f ( x ; μ , σ 2 ) = 1 2 π σ e − ( x − μ ) 2 2 σ 2 fleft(x ; mu, sigma^2right)=frac{1}{sqrt{2 pi} sigma} e^{-frac{(x-mu)^2}{2 sigma^2}} f(x;μ,σ2)=2πσ1e−2σ2(x−μ)2
似然函数为:
L ( μ , σ ∣ x ) = ∏ i = 1 n 1 2 π σ e − ( x i − μ ) 2 2 σ 2 ln L ( μ , σ ∣ x ) = − nln 2 π σ − 1 2 σ 2 ∑ i = 1 n ( x i − μ ) 2 { d d μ ln L ( μ , σ ∣ x ) = 1 2 σ 2 2 ∑ i = 1 n ( x i − μ ) = 0 d d σ ln L ( μ , σ ∣ x ) = − n σ + 1 σ 3 2 ∑ i = 1 n ( x i − μ ) 2 = 0 begin{gathered} L(boldsymbol{mu, sigma} mid boldsymbol{x})=prod_{i=1}^n frac{1}{sqrt{2 pi} sigma} e^{-frac{left(x_i-muright)^2}{2 sigma^2}} \ ln L(boldsymbol{mu, sigma} mid boldsymbol{x})=-operatorname{nln} sqrt{2 pi} sigma-frac{1}{2 sigma^2} sum_{i=1}^nleft(x_i-muright)^2 \ left{begin{array}{c} frac{d}{d mu} ln L(boldsymbol{mu, sigma} mid boldsymbol{x})=frac{1}{2 sigma^2} 2 sum_{i=1}^nleft(x_i-muright)=0 \ frac{d}{d sigma} ln L(boldsymbol{mu, sigma} mid boldsymbol{x})=-frac{n}{sigma}+frac{1}{sigma^3} 2 sum_{i=1}^nleft(x_i-muright)^2=0 end{array}right. end{gathered} L(μ,σ∣x)=i=1∏n2πσ1e−2σ2(xi−μ)2lnL(μ,σ∣x)=−nln2πσ−2σ21i=1∑n(xi−μ)2{dμdlnL(μ,σ∣x)=2σ212∑i=1n(xi−μ)=0dσdlnL(μ,σ∣x)=−σn+σ312∑i=1n(xi−μ)2=0
最后解得:
{ μ = x ˉ σ 2 = 1 n ∑ i = 1 n ( x i − μ ) 2 left{begin{array}{c} mu=bar{x} \ sigma^2=frac{1}{n} sum_{i=1}^nleft(x_i-muright)^2 end{array}right. {μ=xˉσ2=n1∑i=1n(xi−μ)2
因此得到 μ , σ 2 mu, sigma^2 μ,σ2 的最大似然估计量分别为
μ ^ = X ˉ , σ ^ 2 = 1 n ∑ i = 1 n ( X i − X ˉ ) 2 hat{mu}=bar{X}, hat{sigma}^2=frac{1}{n} sum_{i=1}^nleft(X_i-bar{X}right)^2 μ^=Xˉ,σ^2=n1i=1∑n(Xi−Xˉ)2 -
例3 (多项分布): 设 X X X 取值于 { 1 , 2 , 3 , 4 } {1,2,3,4} {1,2,3,4}, 分布概率为 p ( x ∣ θ ) ≜ π x ( θ ) : p(x mid theta) triangleq pi_x(theta): p(x∣θ)≜πx(θ):
π 1 ( θ ) = 1 4 ( 2 + θ ) , π 2 ( θ ) = 1 4 ( 1 − θ ) π 3 ( θ ) = 1 4 ( 1 − θ ) , π 4 ( θ ) = 1 4 θ begin{aligned} pi_1(theta) &=frac{1}{4}(2+theta), & pi_2(theta) &=frac{1}{4}(1-theta) \ pi_3(theta) &=frac{1}{4}(1-theta), & pi_4(theta) &=frac{1}{4} theta end{aligned} π1(θ)π3(θ)=41(2+θ),=41(1−θ),π2(θ)π4(θ)=41(1−θ)=41θ
设 n n n 次试验得到的 X X X 值有 n j n_j nj 个 j ( j = 1 , 2 , 3 , 4 ) j(j=1,2,3,4) j(j=1,2,3,4) ,求参数 θ theta θ 的最大似然估计。解:观测数据的对数似然函数为(去掉了与参数无关的加性常数)
l ( θ ) = n 1 ln ( 2 + θ ) + ( n 2 + n 3 ) ln ( 1 − θ ) + n 4 ln θ . l(theta)=n_1 ln (2+theta)+left(n_2+n_3right) ln (1-theta)+n_4 ln theta . l(θ)=n1ln(2+θ)+(n2+n3)ln(1−θ)+n4lnθ.
令 l ′ ( θ ) = 0 l^{prime}(theta)=0 l′(θ)=0 得到一个关于 θ theta θ 的二次方程,由此写出 argmax l ( θ ) operatorname{argmax} l(theta) argmaxl(θ) 的解析表达式:
l ′ ( θ ) = n 1 2 + θ − n 2 + n 3 1 − θ + n 4 θ l^{prime}(theta)=frac{n_1}{2+theta}-frac{n_2+n_3}{1-theta}+frac{n_4}{theta} l′(θ)=2+θn1−1−θn2+n3+θn4
令 l ′ ( θ ) = 0 l^{prime}(theta)=0 l′(θ)=0 ,即求解二次方程
n θ 2 + ( − n 1 + 2 n 2 + 2 n 3 + n 4 ) θ − 2 n 4 = 0. n theta^2+left(-n_1+2 n_2+2 n_3+n_4right) theta-2 n_4=0 . nθ2+(−n1+2n2+2n3+n4)θ−2n4=0.
得最大似然估计为
θ ^ = − b + b 2 + 8 n n 4 2 n hat{theta}=frac{-b+sqrt{b^2+8 n n_4}}{2 n} θ^=2n−b+b2+8nn4
其中 b = − n 1 + 2 n 2 + 2 n 3 + n 4 b=-n_1+2 n_2+2 n_3+n_4 b=−n1+2n2+2n3+n4 。
2. EM算法
2.1 原理
EM算法最初用于缺失数据模型参数估计,现在已经用在许多优化问题中。设模型中包含
X
obs
和
X
m
i
s
boldsymbol{X}_{text {obs }} 和 boldsymbol{X}_{mathrm{mis}}
Xobs 和Xmis 两个 随机成分, 有联合密度函数或概率函数
f
(
x
o
b
s
,
x
m
i
s
∣
θ
)
,
θ
fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}right), boldsymbol{theta}
f(xobs,xmis∣θ),θ 为未知参数。称
f
(
x
o
b
s
,
x
m
i
s
∣
θ
)
fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}right)
f(xobs,xmis∣θ) 为完全数据的 密度,一般具有简单的形式。实际上我们只有
X
o
b
s
boldsymbol{X}_{mathrm{obs}}
Xobs 的观测数据
X
o
b
s
=
x
o
b
s
,
X
m
i
s
boldsymbol{X}_{mathrm{obs}}=boldsymbol{x}_{mathrm{obs}} , boldsymbol{X}_{mathrm{mis}}
Xobs=xobs,Xmis 不能观测得到, 这一 部分可能是缺失观测数据,也可能是潜在影响因素。所以实际的似然函数为
L
(
θ
)
=
f
(
x
o
b
s
∣
θ
)
=
∫
f
(
x
o
b
s
,
x
m
i
s
∣
θ
)
d
x
m
i
s
,
(1)
tag{1} L(boldsymbol{theta})=fleft(boldsymbol{x}_{mathrm{obs}} mid boldsymbol{theta}right)=int fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}right) d boldsymbol{x}_{mathrm{mis}},
L(θ)=f(xobs∣θ)=∫f(xobs,xmis∣θ)dxmis,(1)
这个似然函数通常比完全数据的似然函数复杂得多,所以很难直接从
L
(
θ
)
L(boldsymbol{theta})
L(θ) 求最大似然估计。
EM算法的想法是,已经有了参数的近似估计值
θ
(
t
)
boldsymbol{theta}^{(t)}
θ(t) 后, 假设
(
X
o
b
s
,
X
m
i
s
)
left(boldsymbol{X}_{mathrm{obs}}, boldsymbol{X}_{mathrm{mis}}right)
(Xobs,Xmis) 近似服从完全密度
f
(
x
o
b
s
,
x
m
i
s
∣
θ
(
t
)
)
fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}^{(t)}right)
f(xobs,xmis∣θ(t)), 这里
X
o
b
s
=
x
o
b
s
boldsymbol{X}_{mathrm{obs}}=boldsymbol{x}_{mathrm{obs}}
Xobs=xobs 已知,所以认为
X
m
i
s
boldsymbol{X}_{mathrm{mis}}
Xmis 近似服从由
f
(
x
o
b
s
,
x
m
i
s
∣
θ
(
t
)
)
fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}^{(t)}right)
f(xobs,xmis∣θ(t)) 导出的条件分 布
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
=
f
(
x
o
b
s
,
x
m
i
s
∣
θ
(
t
)
)
f
(
x
o
b
s
∣
θ
(
t
)
)
(2)
tag{2} fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right)=frac{fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}^{(t)}right)}{fleft(boldsymbol{x}_{mathrm{obs}} mid boldsymbol{theta}^{(t)}right)}
f(xmis∣xobs,θ(t))=f(xobs∣θ(t))f(xobs,xmis∣θ(t))(2)
其中
f
(
x
o
b
s
∣
θ
(
t
)
)
fleft(boldsymbol{x}_{mathrm{obs}} mid boldsymbol{theta}^{(t)}right)
f(xobs∣θ(t)) 是由
f
(
x
o
b
s
,
x
m
i
s
∣
θ
(
t
)
)
fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}^{(t)}right)
f(xobs,xmis∣θ(t)) 决定的边缘密度。据此近似条件分布,在完全数据对数似然函数
ln
f
(
X
o
b
s
,
X
m
i
s
∣
θ
)
ln fleft(boldsymbol{X}_{mathrm{obs}}, boldsymbol{X}_{mathrm{mis}} mid boldsymbol{theta}right)
lnf(Xobs,Xmis∣θ) 中, 把
X
o
b
s
=
x
o
b
s
boldsymbol{X}_{mathrm{obs}}=boldsymbol{x}_{mathrm{obs}}
Xobs=xobs 看成已知, 关于未知部分
X
m
i
s
boldsymbol{X}_{mathrm{mis}}
Xmis 按密度
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{boldsymbol { x }}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right)
f(xmis∣xobs,θ(t)) 求期 望,得到
θ
boldsymbol{theta}
θ 的函数
Q
t
(
θ
)
Q_t(boldsymbol{theta})
Qt(θ) ,再求
Q
t
(
θ
)
Q_t(boldsymbol{theta})
Qt(θ) 的最大值点作为下一个
θ
(
t
+
1
)
boldsymbol{theta}^{(t+1)}
θ(t+1) 。
EM算法每次迭代有如下的E步(期望步)和M步(最大化步):
- E步: 计算完全数据对数似然函数的期望 Q t ( θ ) = E { ln f ( x o b s , X m i s ∣ θ ) } Q_t(boldsymbol{theta})=Eleft{ln fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{X}_{mathrm{mis}} mid boldsymbol{theta}right)right} Qt(θ)=E{lnf(xobs,Xmis∣θ)}, 其中期望针对随机变量 X m i s boldsymbol{X}_{mathrm{mis}} Xmis , 求期望时假定 X m i s boldsymbol{X}_{mathrm{mis}} Xmis 服从条件密度 f ( x m i s ∣ x o b s , θ ( t ) ) fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) f(xmis∣xobs,θ(t)) 决定的分布。
- M步: 求 Q t ( θ ) Q_t(boldsymbol{theta}) Qt(θ) 的最大值点,记为 θ ( t + 1 ) boldsymbol{theta}^{(t+1)} θ(t+1) , 迭代进入下一步。
定理1: EM算法得到的估计序列
θ
(
t
)
boldsymbol{theta}^{(t)}
θ(t) 使得公式(1)中的似然函数值
L
(
θ
(
t
)
)
Lleft(boldsymbol{theta}^{(t)}right)
L(θ(t)) 单调不减。
证明: 对任意参数
θ
boldsymbol{theta}
θ ,有
ln
L
(
θ
)
=
ln
f
(
x
o
b
s
∣
θ
)
⋅
∫
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
=
ln
f
(
x
o
b
s
,
x
m
i
s
∣
θ
)
f
(
x
m
i
s
∣
x
o
b
s
,
θ
)
∫
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
=
∫
[
ln
f
(
x
o
b
s
,
x
m
i
s
∣
θ
)
−
ln
f
(
x
m
i
s
∣
x
o
b
s
,
θ
)
]
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
=
∫
ln
f
(
x
o
b
s
,
x
m
i
s
∣
θ
)
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
−
∫
ln
f
(
x
m
i
s
∣
x
o
b
s
,
θ
)
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
=
Q
t
(
θ
)
−
∫
ln
f
(
x
m
i
s
∣
x
o
b
s
,
θ
)
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
begin{aligned} ln L(theta)=& ln fleft(boldsymbol{x}_{mathrm{obs}} mid thetaright) cdot int fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} \ =&lnfrac{fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}right)}{fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}right)} int fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} \ =& intleft[ln fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}right)-ln fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}right)right] fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} \ =& int ln fleft(boldsymbol{x}_{mathrm{obs}}, boldsymbol{x}_{mathrm{mis}} mid boldsymbol{theta}right) fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} \ &-int ln fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}right) fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} \ =& Q_t(boldsymbol{theta})-int ln fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}right) fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} end{aligned}
lnL(θ)=====lnf(xobs∣θ)⋅∫f(xmis∣xobs,θ(t))dxmislnf(xmis∣xobs,θ)f(xobs,xmis∣θ)∫f(xmis∣xobs,θ(t))dxmis∫[lnf(xobs,xmis∣θ)−lnf(xmis∣xobs,θ)]f(xmis∣xobs,θ(t))dxmis∫lnf(xobs,xmis∣θ)f(xmis∣xobs,θ(t))dxmis−∫lnf(xmis∣xobs,θ)f(xmis∣xobs,θ(t))dxmisQt(θ)−∫lnf(xmis∣xobs,θ)f(xmis∣xobs,θ(t))dxmis
由信息不等式知
∫
ln
f
(
x
m
i
s
∣
x
o
b
s
,
θ
)
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
≤
∫
ln
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
begin{aligned} & int ln fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}right) fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} \ leq & int ln fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} end{aligned}
≤∫lnf(xmis∣xobs,θ)f(xmis∣xobs,θ(t))dxmis∫lnf(xmis∣xobs,θ(t))f(xmis∣xobs,θ(t))dxmis
又EM迭代使得
Q
t
(
θ
(
t
+
1
)
)
≥
Q
t
(
θ
(
t
)
)
Q_tleft(boldsymbol{theta}^{(t+1)}right) geq Q_tleft(boldsymbol{theta}^{(t)}right)
Qt(θ(t+1))≥Qt(θ(t)), 所以
ln
L
(
θ
(
t
+
1
)
)
≥
Q
t
(
θ
(
t
)
)
−
∫
ln
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
f
(
x
m
i
s
∣
x
o
b
s
,
θ
(
t
)
)
d
x
m
i
s
=
ln
L
(
θ
(
t
)
)
.
begin{aligned} ln Lleft(boldsymbol{theta}^{(t+1)}right) & geq Q_tleft(boldsymbol{theta}^{(t)}right)-int ln fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) fleft(boldsymbol{x}_{mathrm{mis}} mid boldsymbol{x}_{mathrm{obs}}, boldsymbol{theta}^{(t)}right) d boldsymbol{x}_{mathrm{mis}} \ &=ln Lleft(boldsymbol{theta}^{(t)}right) . end{aligned}
lnL(θ(t+1))≥Qt(θ(t))−∫lnf(xmis∣xobs,θ(t))f(xmis∣xobs,θ(t))dxmis=lnL(θ(t)).
定理证毕。
在适当正则性条件下, EM算法的迭代序列依概率收敛到 θ ( t ) theta^{(t)} θ(t)的最大值点。 但是, 定理(1)仅保证EM算法最终能收敛, 但不能保证EM算法会收敛到似然函数的全局最大值点, 算法也可能收敛到局部极大值点或者鞍点。
在实际问题中, 往往E步和M步都比较简单, 有时E步和M步都有解析表达式, 这时EM算法实现很简单。 EM算法优点是计算稳定, 可以保持原有的参数约束, 缺点是收敛可能很慢, 尤其是接近最大值点时可能收敛更慢。 如果公式(1)中的似然函数不是凸函数, 算法可能收敛不到全局最大值点, 遇到这样的问题可以多取不同初值比较, 用矩估计等合适的近似值作为初值。
原理暂时看不懂没关系,结合后面例题看就更容易懂了。
2.2 示例
(混合分布) EM算法可以用来估计混合分布的参数。 设随机变量
Y
1
∼
N
(
μ
1
,
δ
1
)
,
Y
2
∼
N
(
μ
2
,
δ
2
)
Y_1 sim mathrm{N}left(mu_1, delta_1right), Y_2 sim mathrm{N}left(mu_2, delta_2right)
Y1∼N(μ1,δ1),Y2∼N(μ2,δ2),
Y
1
,
Y
2
Y_1, Y_2
Y1,Y2 独立。记
N
(
μ
,
δ
)
N(mu, delta)
N(μ,δ) 的密度为
f
(
x
∣
μ
,
δ
)
f(x mid mu, delta)
f(x∣μ,δ) 。设随机变量
W
∼
b
(
1
,
λ
)
,
0
<
λ
<
1
,
W
W sim mathrm{b}(1, lambda), 0<lambda<1, W
W∼b(1,λ),0<λ<1,W 与
Y
1
,
Y
2
Y_1, Y_2
Y1,Y2 独立,令
X
=
(
1
−
W
)
Y
1
+
W
Y
2
,
X=(1-W) Y_1+W Y_2,
X=(1−W)Y1+WY2,
则
W
=
0
W=0
W=0 条件下
X
∼
N
(
μ
1
,
δ
1
)
,
W
=
1
X sim mathrm{N}left(mu_1, delta_1right), W=1
X∼N(μ1,δ1),W=1 条件下
X
∼
N
(
μ
2
,
δ
2
)
X sim mathrm{N}left(mu_2, delta_2right)
X∼N(μ2,δ2), 但
X
X
X 的边缘密度为
f
(
x
∣
θ
)
=
(
1
−
λ
)
f
(
x
∣
μ
1
,
δ
1
)
+
λ
f
(
x
∣
μ
2
,
δ
2
)
,
f(x mid boldsymbol{theta})=(1-lambda) fleft(x mid mu_1, delta_1right)+lambda fleft(x mid mu_2, delta_2right),
f(x∣θ)=(1−λ)f(x∣μ1,δ1)+λf(x∣μ2,δ2),
其中
θ
=
(
μ
1
,
δ
1
,
μ
2
,
δ
2
,
λ
)
boldsymbol{theta}=left(mu_1, delta_1, mu_2, delta_2, lambdaright)
θ=(μ1,δ1,μ2,δ2,λ) 。
设
X
X
X 有样本
X
=
(
X
1
,
…
,
X
n
)
boldsymbol{X}=left(X_1, ldots, X_nright)
X=(X1,…,Xn), 样本值为
x
boldsymbol{x}
x , 实际观测数据的似然函数为
L
(
θ
)
=
∏
i
=
1
n
f
(
x
i
∣
θ
)
L(boldsymbol{theta})=prod_{i=1}^n fleft(x_i mid boldsymbol{theta}right)
L(θ)=i=1∏nf(xi∣θ)
这个函数是光滑函数但是形状很复杂, 直接求极值很容易停留在局部极值点。
用EM算法,以
W
=
(
W
1
,
…
,
W
n
)
boldsymbol{W}=left(W_1, ldots, W_nright)
W=(W1,…,Wn) 为没有观测到的部分, 完全数据的似然函数和对数似然函数为
L
~
(
θ
∣
x
,
W
)
=
∏
W
i
=
0
f
(
x
i
∣
μ
1
,
δ
1
)
∏
W
i
=
1
f
(
x
i
∣
μ
2
,
δ
2
)
λ
∑
i
=
1
n
W
i
(
1
−
λ
)
n
−
∑
i
=
1
n
W
i
l
~
(
θ
∣
x
,
W
)
=
∑
i
=
1
n
[
(
1
−
W
i
)
log
f
(
x
i
∣
μ
1
,
δ
1
)
+
W
i
log
f
(
x
i
∣
μ
2
,
δ
2
)
]
+
(
∑
i
=
1
n
W
i
)
log
λ
+
(
n
−
∑
i
=
1
n
W
i
)
log
(
1
−
λ
)
begin{aligned} tilde{L}(boldsymbol{theta} mid boldsymbol{x}, boldsymbol{W})=& prod_{W_i=0} fleft(x_i mid mu_1, delta_1right) prod_{W_i=1} fleft(x_i mid mu_2, delta_2right) lambda^{sum_{i=1}^n W_i}(1-lambda)^{n-sum_{i=1}^n W_i} \ tilde{l}(boldsymbol{theta} mid boldsymbol{x}, boldsymbol{W})=& sum_{i=1}^nleft[left(1-W_iright) log fleft(x_i mid mu_1, delta_1right)+W_i log fleft(x_i mid mu_2, delta_2right)right] \ &+left(sum_{i=1}^n W_iright) log lambda+left(n-sum_{i=1}^n W_iright) log (1-lambda) end{aligned}
L~(θ∣x,W)=l~(θ∣x,W)=Wi=0∏f(xi∣μ1,δ1)Wi=1∏f(xi∣μ2,δ2)λ∑i=1nWi(1−λ)n−∑i=1nWii=1∑n[(1−Wi)logf(xi∣μ1,δ1)+Wilogf(xi∣μ2,δ2)]+(i=1∑nWi)logλ+(n−i=1∑nWi)log(1−λ)
在E步,设已有
θ
boldsymbol{theta}
θ 的近似值
θ
(
t
)
=
(
μ
1
(
t
)
,
δ
1
(
t
)
,
μ
2
(
t
)
,
δ
2
(
t
)
,
λ
(
t
)
)
boldsymbol{theta}^{(t)}=left(mu_1^{(t)}, delta_1^{(t)}, mu_2^{(t)}, delta_2^{(t)}, lambda^{(t)}right)
θ(t)=(μ1(t),δ1(t),μ2(t),δ2(t),λ(t)), 以
θ
(
t
)
boldsymbol{theta}^{(t)}
θ(t) 为分布参数,在
X
=
x
boldsymbol{X}=boldsymbol{x}
X=x 条件下,
W
i
boldsymbol{W}_i
Wi 的 条件分布为
γ
i
(
t
)
≜
P
(
W
i
=
1
∣
x
,
θ
(
t
)
)
=
P
(
W
i
=
1
∣
X
i
=
x
i
,
θ
(
t
)
)
=
λ
(
t
)
f
(
x
i
∣
μ
2
(
t
)
,
δ
2
(
t
)
)
(
1
−
λ
(
t
)
)
f
(
x
i
∣
μ
1
(
t
)
,
δ
1
(
t
)
)
+
λ
(
t
)
f
(
x
i
∣
μ
2
(
t
)
,
δ
2
(
t
)
)
.
begin{aligned} gamma_i^{(t)} & triangleq Pleft(W_i=1 mid boldsymbol{x}, boldsymbol{theta}^{(t)}right)=Pleft(W_i=1 mid X_i=x_i, boldsymbol{theta}^{(t)}right) \ &=frac{lambda^{(t)} fleft(x_i mid mu_2^{(t)}, delta_2^{(t)}right)}{left(1-lambda^{(t)}right) fleft(x_i mid mu_1^{(t)}, delta_1^{(t)}right)+lambda^{(t)} fleft(x_i mid mu_2^{(t)}, delta_2^{(t)}right)} . end{aligned}
γi(t)≜P(Wi=1∣x,θ(t))=P(Wi=1∣Xi=xi,θ(t))=(1−λ(t))f(xi∣μ1(t),δ1(t))+λ(t)f(xi∣μ2(t),δ2(t))λ(t)f(xi∣μ2(t),δ2(t)).
这里的推导类似于逆概率公式。利用
W
i
W_i
Wi 的条件分布求完全数据对数似然的期望,得
Q
t
(
θ
)
=
∑
i
=
1
n
[
(
1
−
γ
i
(
t
)
)
log
f
(
x
i
∣
μ
1
,
δ
1
)
+
γ
i
(
t
)
log
f
(
x
i
∣
μ
2
,
δ
2
)
]
+
(
∑
i
=
1
n
γ
i
(
t
)
)
log
λ
+
(
n
−
∑
i
=
1
n
γ
i
(
t
)
)
log
(
1
−
λ
)
.
(3)
tag{3} begin{aligned} Q_t(boldsymbol{theta})=& sum_{i=1}^nleft[left(1-gamma_i^{(t)}right) log fleft(x_i mid mu_1, delta_1right)+gamma_i^{(t)} log fleft(x_i mid mu_2, delta_2right)right] \ &+left(sum_{i=1}^n gamma_i^{(t)}right) log lambda+left(n-sum_{i=1}^n gamma_i^{(t)}right) log (1-lambda) . end{aligned}
Qt(θ)=i=1∑n[(1−γi(t))logf(xi∣μ1,δ1)+γi(t)logf(xi∣μ2,δ2)]+(i=1∑nγi(t))logλ+(n−i=1∑nγi(t))log(1−λ).(3)
令
∇
Q
t
(
θ
)
=
0
nabla Q_t(boldsymbol{theta})=mathbf{0}
∇Qt(θ)=0 ,求得
Q
t
(
θ
)
Q_t(boldsymbol{theta})
Qt(θ) 的最大值点
θ
(
t
+
1
)
boldsymbol{theta}^{(t+1)}
θ(t+1) 为
{
μ
1
(
t
+
1
)
=
∑
i
=
1
n
(
1
−
γ
i
(
t
)
)
x
i
∑
i
=
1
n
(
1
−
γ
i
(
t
)
)
δ
1
(
t
+
1
)
=
∑
i
=
1
n
(
1
−
γ
i
(
t
)
)
(
x
i
−
μ
1
(
t
+
1
)
)
2
∑
i
=
1
n
(
1
−
γ
i
(
t
)
)
μ
2
(
t
+
1
)
=
∑
i
=
1
n
γ
i
(
t
)
x
i
∑
i
=
1
n
γ
i
(
t
)
δ
2
(
t
+
1
)
=
∑
i
=
1
n
γ
i
(
t
)
(
x
i
−
μ
2
(
t
+
1
)
)
2
∑
i
=
1
n
γ
i
(
t
)
λ
(
t
+
1
)
=
1
n
∑
i
=
1
n
γ
i
(
t
)
(4)
tag{4} left{begin{array}{l} mu_1^{(t+1)}=frac{sum_{i=1}^nleft(1-gamma_i^{(t)}right) x_i}{sum_{i=1}^nleft(1-gamma_i^{(t)}right)} \ delta_1^{(t+1)}=frac{sum_{i=1}^nleft(1-gamma_i^{(t)}right)left(x_i-mu_1^{(t+1)}right)^2}{sum_{i=1}^nleft(1-gamma_i^{(t)}right)} \ mu_2^{(t+1)}=frac{sum_{i=1}^n gamma_i^{(t)} x_i}{sum_{i=1}^n gamma_i^{(t)}} \ delta_2^{(t+1)}=frac{sum_{i=1}^n gamma_i^{(t)}left(x_i-mu_2^{(t+1)}right)^2}{sum_{i=1}^n gamma_i^{(t)}} \ lambda^{(t+1)}=frac{1}{n} sum_{i=1}^n gamma_i^{(t)} end{array}right.
⎩
⎨
⎧μ1(t+1)=∑i=1n(1−γi(t))∑i=1n(1−γi(t))xiδ1(t+1)=∑i=1n(1−γi(t))∑i=1n(1−γi(t))(xi−μ1(t+1))2μ2(t+1)=∑i=1nγi(t)∑i=1nγi(t)xiδ2(t+1)=∑i=1nγi(t)∑i=1nγi(t)(xi−μ2(t+1))2λ(t+1)=n1∑i=1nγi(t)(4)
适当选取初值
θ
(
0
)
boldsymbol{theta}^{(0)}
θ(0) 用公式(3)和(4)迭代就可以计算
θ
boldsymbol{theta}
θ 的最大似然估计。
-
(多项分布): 设 X X X 取值于 { 1 , 2 , 3 , 4 } {1,2,3,4} {1,2,3,4}, 分布概率为 p ( x ∣ θ ) ≜ π x ( θ ) : p(x mid theta) triangleq pi_x(theta): p(x∣θ)≜πx(θ):
π 1 ( θ ) = 1 4 ( 2 + θ ) , π 2 ( θ ) = 1 4 ( 1 − θ ) π 3 ( θ ) = 1 4 ( 1 − θ ) , π 4 ( θ ) = 1 4 θ begin{aligned} pi_1(theta) &=frac{1}{4}(2+theta), & pi_2(theta) &=frac{1}{4}(1-theta) \ pi_3(theta) &=frac{1}{4}(1-theta), & pi_4(theta) &=frac{1}{4} theta end{aligned} π1(θ)π3(θ)=41(2+θ),=41(1−θ),π2(θ)π4(θ)=41(1−θ)=41θ
设 n n n 次试验得到的 X X X 值有 n j n_j nj 个 j ( j = 1 , 2 , 3 , 4 ) j(j=1,2,3,4) j(j=1,2,3,4) ,求参数 θ theta θ 的最大似然估计。P ( X = 12 ) = 1 4 θ 。 P(X=12)=frac{1}{4} theta_{text {。 }} P(X=12)=41θ。 这时 Z 2 + n 4 Z_2+n_4 Z2+n4 代表结果12和结果4的出现次数, 这两种结果出现概率为 1 2 θ frac{1}{2} theta 21θ ,其它结 果 ( 11 , 2 , 3 ) (11,2,3) (11,2,3) 的出现概率为 1 − 1 2 θ 。 1-frac{1}{2} theta_{text {。 }} 1−21θ。 令 Y = Z 2 + n 4 Y=Z_2+n_4 Y=Z2+n4 ,则 Y ∼ B ( n , 1 2 θ ) Y sim mathrm{B}left(n, frac{1}{2} thetaright) Y∼B(n,21θ) 。
数据 ( Z 1 , Z 2 , n 2 , n 3 , n 4 ) left(Z_1, Z_2, n_2, n_3, n_4right) (Z1,Z2,n2,n3,n4) 的全似然函数为
L c ( θ ) ∝ ( 1 2 ) Z 1 ( θ 4 ) Z 2 ( 1 4 − θ 4 ) n 2 + n 3 ( θ 4 ) n 4 L_c(theta) proptoleft(frac{1}{2}right)^{Z_1}left(frac{theta}{4}right)^{Z_2}left(frac{1}{4}-frac{theta}{4}right)^{n_2+n_3}left(frac{theta}{4}right)^{n_4} Lc(θ)∝(21)Z1(4θ)Z2(41−4θ)n2+n3(4θ)n4
对数似然函数 (差一个与 θ theta θ 无关的常数项) 为
l c ( θ ) = ln L c ( θ ) = ( Z 2 + n 4 ) ln θ + ( n 2 + n 3 ) ln ( 1 − θ ) l_c(theta)=ln L_c(theta)=left(Z_2+n_4right) ln theta+left(n_2+n_3right) ln (1-theta) lc(θ)=lnLc(θ)=(Z2+n4)lnθ+(n2+n3)ln(1−θ)
在EM迭代中, 假设已经得到的参数 θ theta θ 近似值为 θ ( t ) theta^{(t)} θ(t) , 设 θ = θ ( t ) theta=theta^{(t)} θ=θ(t) , 在给定 n , n 1 , n 2 , n 3 , n 4 n, n_1, n_2, n_3, n_4 n,n1,n2,n3,n4 条件下求 l c ( θ ) l_c(theta) lc(θ) 的 条件期望, 这时 Z 2 Z_2 Z2 的条件分布为
B ( n 1 , θ ( t ) 2 + θ ( t ) ) , mathrm{B}left(n_1, frac{theta^{(t)}}{2+theta^{(t)}}right), B(n1,2+θ(t)θ(t)),
于是
Z 2 ( t ) = E ( Z 2 ∣ θ ( t ) , n , n 1 , n 2 , n 3 , n 4 ) = n 1 θ ( t ) 2 + θ ( t ) , Z_2^{(t)}=Eleft(Z_2 mid theta^{(t)}, n, n_1, n_2, n_3, n_4right)=n_1 frac{theta^{(t)}}{2+theta^{(t)}}, Z2(t)=E(Z2∣θ(t),n,n1,n2,n3,n4)=n12+θ(t)θ(t),
从而完全对数似然函数的条件期望为
Q t ( θ ) = E ( ln L c ( θ ) ∣ θ ( t ) , n , n 1 , n 2 , n 3 , n 4 ) = ( Z 2 ( t ) + n 4 ) ln θ + ( n 2 + n 3 ) ln ( 1 − θ ) . Q_t(theta)=Eleft(ln L_c(theta) mid theta^{(t)}, n, n_1, n_2, n_3, n_4right)=left(Z_2^{(t)}+n_4right) ln theta+left(n_2+n_3right) ln (1-theta) . Qt(θ)=E(lnLc(θ)∣θ(t),n,n1,n2,n3,n4)=(Z2(t)+n4)lnθ+(n2+n3)ln(1−θ).
求解 Q t ( θ ) Q_t(theta) Qt(θ) 的最大值,令
d d θ Q t ( θ ) = n 4 + Z 2 ( t ) θ − n 2 + n 3 1 − θ = 0 frac{d}{d theta} Q_t(theta)=frac{n_4+Z_2^{(t)}}{theta}-frac{n_2+n_3}{1-theta}=0 dθdQt(θ)=θn4+Z2(t)−1−θn2+n3=0
得下一个参数近似值为
θ ( t + 1 ) = n 4 + Z 2 ( t ) n 2 + n 3 + n 4 + Z 2 ( t ) . theta^{(t+1)}=frac{n_4+Z_2^{(t)}}{n_2+n_3+n_4+Z_2^{(t)}} . θ(t+1)=n2+n3+n4+Z2(t)n4+Z2(t).
于是, EM迭代步骤从某个 θ ( 0 ) theta^{(0)} θ(0) 出发,比如 θ ( 0 ) = 1 2 theta^{(0)}=frac{1}{2} θ(0)=21 , 在第 t t t 步计算
Z 2 ( t ) = n 1 θ ( t ) 2 + θ ( t ) , θ ( t + 1 ) = n 4 + Z 2 ( t ) n 2 + n 3 + n 4 + Z 2 ( t ) Z_2^{(t)}=n_1 frac{theta^{(t)}}{2+theta^{(t)}}, quad theta^{(t+1)}=frac{n_4+Z_2^{(t)}}{n_2+n_3+n_4+Z_2^{(t)}} Z2(t)=n12+θ(t)θ(t),θ(t+1)=n2+n3+n4+Z2(t)n4+Z2(t)
迭代到两次的近似参数值变化小于 ϵ = 1 0 − 6 epsilon=10^{-6} ϵ=10−6 为止。
最后
以上就是自然蓝天为你收集整理的统计学 | 最大似然估计与EM算法(持续更新)参考资料1. 最大似然估计2. EM算法的全部内容,希望文章能够帮你解决统计学 | 最大似然估计与EM算法(持续更新)参考资料1. 最大似然估计2. EM算法所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复