我是靠谱客的博主 自然蓝天,最近开发中收集的这篇文章主要介绍统计学 | 最大似然估计与EM算法(持续更新)参考资料1. 最大似然估计2. EM算法,觉得挺不错的,现在分享给大家,希望可以做个参考。

概述

文章目录

  • 参考资料
  • 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=1np(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=1nlnp(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 Xb(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(1p)1x
    故似然函数为:
    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(px)=i=1npxi(1p)1xi=pi=1nxi(1p)ni=1nxilnL(px)=i=1nxilnp+(ni=1nxi)ln(1p)
    对待估参数求导为 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(px)=pi=1nxi1pni=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 XN(μ,σ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π σ1e2σ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=1n2π σ1e2σ2(xiμ)2lnL(μ,σx)=nln2π σ2σ21i=1n(xiμ)2{dμdlnL(μ,σx)=2σ212i=1n(xiμ)=0dσdlnL(μ,σx)=σn+σ312i=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=n1i=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=1n(XiXˉ)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+θn11θ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} θ^=2nb+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=xobsXmis 不能观测得到, 这一 部分可能是缺失观测数据,也可能是潜在影响因素。所以实际的似然函数为
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(xmisxobs,θ(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(xmisxobs,θ(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(xmisxobs,θ(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(xmisxobs,θ(t))dxmislnf(xmisxobs,θ)f(xobs,xmisθ)f(xmisxobs,θ(t))dxmis[lnf(xobs,xmisθ)lnf(xmisxobs,θ)]f(xmisxobs,θ(t))dxmislnf(xobs,xmisθ)f(xmisxobs,θ(t))dxmislnf(xmisxobs,θ)f(xmisxobs,θ(t))dxmisQt(θ)lnf(xmisxobs,θ)f(xmisxobs,θ(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(xmisxobs,θ)f(xmisxobs,θ(t))dxmislnf(xmisxobs,θ(t))f(xmisxobs,θ(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(xmisxobs,θ(t))f(xmisxobs,θ(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) Y1N(μ1,δ1),Y2N(μ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 Wb(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=(1W)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 XN(μ1,δ1),W=1 条件下 X ∼ N ( μ 2 , δ 2 ) X sim mathrm{N}left(mu_2, delta_2right) XN(μ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=1nf(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=0f(xiμ1,δ1)Wi=1f(xiμ2,δ2)λi=1nWi(1λ)ni=1nWii=1n[(1Wi)logf(xiμ1,δ1)+Wilogf(xiμ2,δ2)]+(i=1nWi)logλ+(ni=1nWi)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=1x,θ(t))=P(Wi=1Xi=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=1n[(1γi(t))logf(xiμ1,δ1)+γi(t)logf(xiμ2,δ2)]+(i=1nγi(t))logλ+(ni=1nγ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)=n1i=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 {。 }} 121θ  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) YB(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(414θ)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} ϵ=106 为止。

最后

以上就是自然蓝天为你收集整理的统计学 | 最大似然估计与EM算法(持续更新)参考资料1. 最大似然估计2. EM算法的全部内容,希望文章能够帮你解决统计学 | 最大似然估计与EM算法(持续更新)参考资料1. 最大似然估计2. EM算法所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。

本图文内容来源于网友提供,作为学习参考使用,或来自网络收集整理,版权属于原作者所有。
点赞(49)

评论列表共有 0 条评论

立即
投稿
返回
顶部