谈谈期望最大算法

2019年2月13日12:28:39 发表评论浏览:32

本篇讨论期望最大算法算法(Expecation Maximization Algorithm,EM 算法)。

从极大似然估计到 EM 算法

EM 算法和极大似然有着很深的联系。因此在讨论 EM 算法之前,我们首先通过一个例子来回顾极大似然算法。

极大似然

问题描述

假设有一枚硬币 A;投掷硬币 A 后结果为正面的概率记为 pp。我们将投掷硬币结果为正面的实验记作 11,反面的实验记作 00。独立重复 n=10n=10 次实验,观测结果如下:1,1,0,1,0,0,1,0,1,1.1,1,0,1,0,0,1,0,1,1.

现在的问题是要估计投掷硬币 A 后结果为正面的概率 pp。

整理似然函数和对数似然函数

显而易见,投掷硬币 A 后的结果应当服从二项分布。于是我们可以整理似然函数:L(θ)=L(x1,x2,…,x10∣θ)=∏i=110p(xi∣θ)=p6(1−p)4.L(θ)=L(x1,x2,…,x10∣θ)=∏i=110p(xi∣θ)=p6(1−p)4.

由于连乘不方便分析和优化,因此我们对似然函数取对数,构造对数似然函数,将连乘转换为连加。H(θ)=lnL(θ)=∑i=110lnp(xi∣θ)=(lnp)6⋅(ln(1−p))4.H(θ)=ln⁡L(θ)=∑i=110ln⁡p(xi∣θ)=(ln⁡p)6⋅(ln⁡(1−p))4.

求解极大化问题

根据似然函数或者对数似然函数,构造极大化问题。此处以似然函数为例:θ^=argmaxθL(θ).θ^=argmaxθL(θ).

对于这类简单的似然函数,我们可以对它求导,令导数为 00,而后得到参数值;对于复杂的似然函数极大化问题,我们可以用牛顿法、拟牛顿法、梯度上升法等方法求解。此处dL(θ)dθ=dLdp=2p5(1−p)3(3−5p).dL(θ)dθ=dLdp=2p5(1−p)3(3−5p).

令导函数为 00,得到似然方程2p5(1−p)3(3−5p)=02p5(1−p)3(3−5p)=0

对似然方程的三个不重复解代入原函数依次验证,可知当 p=35p=35 时原函数取得极大值。此即为所求。

总结

极大似然估计的一般方法可以记录如下:

  1. 极大似然估计,估计的目标是随机变量服从分布的参数;因此在使用极大似然估计法之前,必须假定数据服从的分布。
  2. 根据分布和实验结果,列出似然函数或者对数似然函数;
  3. 求解似然函数或者对数似然函数的极大化问题;
  4. 对于不同的(对数)似然函数,采取不同的求解方法——例如可以采用牛顿法、拟牛顿法、梯度上升法等方法,也可以求导直接求解。

EM 算法

问题描述

我们考虑一个复杂一些的问题。

现在假设我们有 A, B, C 三枚硬币;三枚硬币投掷后结果为正面的概率分别记为 ππ, pp, qq。现在我们规定这样一种实验:先抛 A 硬币,若 A 硬币为正面,则继续抛 B 硬币并记录结果;若 A 硬币为反面,则继续抛 C 硬币并记录结果。独立重复 n=10n=10 次实验,观测结果如下:1,1,0,1,0,0,1,0,1,1.1,1,0,1,0,0,1,0,1,1.

现在的问题是要估计三枚硬币投掷后结果为正面的概率 ππ, pp, qq。

问题分析、生成模型与似然函数

参考上一节关于极大似然估计法的讨论。如果我们能知道 10 次独立重复实验中,哪一些是由 B 硬币确定的最终结果,哪一些是由 C 硬币确定的最终结果;那么我们就可以对他们分别采用极大似然估计,确定 B、C 两枚硬币的二项分布参数。但是由于硬币 A 的分布参数不确定,反过来我们又必须先知道 B、C 两枚硬币取得正面的分布参数,才能说某一个结果更可能来自 B 硬币还是 C 硬币。

这是一个鸡生蛋蛋生鸡的问题,看起来不好解决。脑壳疼啊脑壳疼。以下我们用概率的语言来表述这个问题。

我们引入两个随机变量 yy 和 zz。此处 yy 表示一次投币实验的结果,而 zz 表示一次投币实验中间硬币 A 的结果。于是,三硬币问题的生成模型是p(y∣θ)===∑zp(y,z∣θ)∑zp(z∣θ)p(y∣z,θ)πpy(1−p)1−y+(1−π)qy(1−q)1−y.(*)(*)p(y∣θ)=∑zp(y,z∣θ)=∑zp(z∣θ)p(y∣z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y.

如此有似然函数L(θ)==p(Y∣θ)∏i=1nπpyi(1−p)1−yi+(1−π)qyi(1−q)1−yi.L(θ)=p(Y∣θ)=∏i=1nπpyi(1−p)1−yi+(1−π)qyi(1−q)1−yi.

考虑求模型参数 θ=(π,p,q)θ=(π,p,q) 的极大似然估计θ^=argmaxθlnL(θ).θ^=argmaxθln⁡L(θ).

这种「鸡生蛋蛋生鸡」的问题,(π,p,q)(π,p,q) 参数之间互相制约,因此这个优化问题很难求得解析解。因此,针对极大似然估计法的常用解法在这里都会失效,而只能用迭代的方式求解。EM 算法是其中一种迭代算法。

EM 迭代

我们刚才说,如果能确定三硬币模型的参数,那么我们就能知道观测数据 yiyi 来自硬币 B 或 C 的概率。知道了这个概率,我们又可以反过来去推算三硬币模型的参数。这种鸡生蛋蛋生鸡的问题,总得有个头啊。因此,在 EM 算法迭代过程中,我们需要为模型参数设置一个初始值,而后不断循环迭代。以下是以三硬币模型为例,表述 EM 算法的迭代过程。

  1. 选取参数的初始值,记为 θ(0)=(π(0),p(0),q(0))θ(0)=(π(0),p(0),q(0));
  2. 迭代至参数值变化不大为止
    1. E 步,根据第 i−1i−1 轮的参数,计算隐变量 zz 的概率 μ(i)j=p(zj∣θ(i−1))μj(i)=p(zj∣θ(i−1))(即 yjyj 来自硬币 B 的概率);
    2. M 步,根据隐变量的概率,计算模型参数新的估计值π(i)=p(i)=q(i)=1n∑j=1nμ(i)j;∑nj=1μ(i)jyj∑nj=1μ(i)j;∑nj=1(1−μ(i)j)yj∑nj=1(1−μ(i)j).π(i)=1n∑j=1nμj(i);p(i)=∑j=1nμj(i)yj∑j=1nμj(i);q(i)=∑j=1n(1−μj(i))yj∑j=1n(1−μj(i)).

实际计算看看

我们假定初值 θ(0)=(π(0)=0.5,p(0)=0.5,q(0)=0.5)θ(0)=(π(0)=0.5,p(0)=0.5,q(0)=0.5)。

第一轮迭代:

  • E 步:μ(1)j=p(zj∣θ(0))=π(0)(p(0))yj(1−p(0))1−yjπ(0)(p(0))yj(1−p(0))1−yj+(1−π(0))(q(0))yj(1−q(0))1−yj=0.5μ(1)j=p(zj∣θ(0))=π(0)(p(0))yj(1−p(0))1−yjπ(0)(p(0))yj(1−p(0))1−yj+(1−π(0))(q(0))yj(1−q(0))1−yj=0.5
  • M 步:π(1)=p(1)=q(1)=1n∑j=1nμ(1)j=0.5;∑nj=1μ(1)jyj∑nj=1μ(1)j=0.6;∑nj=1(1−μ(1)j)yj∑nj=1(1−μ(1)j)=0.6.π(1)=1n∑j=1nμj(1)=0.5;p(1)=∑j=1nμj(1)yj∑j=1nμj(1)=0.6;q(1)=∑j=1n(1−μj(1))yj∑j=1n(1−μj(1))=0.6.

第二轮迭代:

  • E 步:μ(2)j=p(zj∣θ(1))=π(1)(p(1))yj(1−p(1))1−yjπ(1)(p(1))yj(1−p(1))1−yj+(1−π(1))(q(1))yj(1−q(1))1−yj=0.5μ(2)j=p(zj∣θ(1))=π(1)(p(1))yj(1−p(1))1−yjπ(1)(p(1))yj(1−p(1))1−yj+(1−π(1))(q(1))yj(1−q(1))1−yj=0.5
  • M 步:π(2)=p(2)=q(2)=1n∑j=1nμ(2)j=0.5;∑nj=1μ(2)jyj∑nj=1μ(2)j=0.6;∑nj=1(1−μ(2)j)yj∑nj=1(1−μ(2)j)=0.6.π(2)=1n∑j=1nμj(2)=0.5;p(2)=∑j=1nμj(2)yj∑j=1nμj(2)=0.6;q(2)=∑j=1n(1−μj(2))yj∑j=1n(1−μj(2))=0.6.

考虑到 θ(1)=θ(2)=(π=0.5,p=0.6,q=0.6)θ(1)=θ(2)=(π=0.5,p=0.6,q=0.6),参数已趋于稳定,故而迭代结束。

注意,EM 算法迭代的结果是对初值敏感的。例如,对于三硬币模型,若初值取 θ(0)=(π(0)=0.4,p(0)=0.6,q(0)=0.7)θ(0)=(π(0)=0.4,p(0)=0.6,q(0)=0.7),则迭代结果为 θ^=(π^=0.4064,p^=0.5368,q^=0.6432)θ^=(π^=0.4064,p^=0.5368,q^=0.6432)。

EM 算法的理论推导

数学工具

凸函数

定义在实数域上的一元函数 f(x)f(x),若对于定义域内的任意自变量都成立不等式 f′′(x)⩾0f″(x)⩾0,则称该函数是凸函数;若严格成立不等式 f′′(x)>0f″(x)>0,则称该函数是严格凸函数。

定义在实数域上的多元函数 f(x1,x2,…,xn)f(x1,x2,…,xn),若对于定义域内的任意自变量其海森矩阵是半正定的,则称该函数是凸函数;若海森矩阵是正定的,则称该函数是严格凸函数。

期望

对于连续随机变量 XX,若其概率密度函数为 f(x)f(x),则其期望定义为E[X]=∫+∞−∞x⋅f(x)dx.E[X]=∫−∞+∞x⋅f(x)dx.

设随机变量 Y=g(X)Y=g(X),则其期望是E[Y]=∫+∞−∞g(x)⋅f(x)dx.E[Y]=∫−∞+∞g(x)⋅f(x)dx.

琴生(Jensen)不等式

假设 μμ 是集合 ΩΩ 的正测度,并且 μ(Ω)=1μ(Ω)=1。又设 gg 是勒贝格(Lebesgue)可积的实值函数,而 φφ 是定义在 RR 上的凸函数,则成立琴生不等式φ(∫Ωgdμ)⩽∫Ωφ∘gdμ.φ(∫Ωgdμ)⩽∫Ωφ∘gdμ.

以概率论的语言来说,μμ 显然是个概率测度。此时,若 gg 表示随机变量 XX,则在 ΩΩ 上,随机变量相对概率测度的积分其实是期望。此时,琴生不等式转化为φ(E[X])⩽E(φ[X]).φ(E[X])⩽E(φ[X]).

问题的数学表述

对于 mm 个相互独立的样本 x1x1, x2x2, …, xmxm,有观测数据 y1y1, y2y2, …, ymym,以及隐含数据 z1z1, z2z2, …, zmzm。此时,(y)(y) 为不完全数据,而 (y,z)(y,z) 是完全数据。

面对不完全数据建立的概率模型(** 式)会含有隐变量:H(θ)====L(θ)lnp(Y∣θ)ln∑zp(Y,Z∣θ)ln(∑zp(Z∣θ)p(Y∣Z,θ))H(θ)=L(θ)=ln⁡p(Y∣θ)=ln⁡∑zp(Y,Z∣θ)=ln⁡(∑zp(Z∣θ)p(Y∣Z,θ))

因此,极大化问题事实上变成了θ^,z^=argmaxθ,zH(θ,z).θ^,z^=argmaxθ,zH(θ,z).

此时,若以极大似然估计的解法思路,对 θθ 和 zz 求导,而后直接求解极大似然方程或者梯度上升,其表达式会非常复杂。因此,很难求得它的解析解。

构造似然函数下界

我们的目的是极大化(对数)似然函数 H(θ)H(θ)。因此对于第 i+1i+1 轮迭代,我们希望找到新的 θθ 使得 H(θ)>H(θ(i))H(θ)>H(θ(i))。

因此,我们考虑(凹函数的)琴生不等式H(θ)==⩾ln(∑zp(Z∣θ)p(Y∣Z,θ))ln(∑zp(Z∣Y,θ(i))⋅p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i)))∑zp(Z∣Y,θ(i))ln(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))).H(θ)=ln⁡(∑zp(Z∣θ)p(Y∣Z,θ))=ln⁡(∑zp(Z∣Y,θ(i))⋅p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i)))⩾∑zp(Z∣Y,θ(i))ln⁡(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))).

注意到 ∑zp(Z∣Y,θ(i))=1∑zp(Z∣Y,θ(i))=1,于是考虑做差H(θ)−H(θ(i))⩾=∑zp(Z∣Y,θ(i))ln(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i)))−lnp(Y∣θ(i))∑zp(Z∣Y,θ(i))ln(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))p(Y∣θ(i))).H(θ)−H(θ(i))⩾∑zp(Z∣Y,θ(i))ln⁡(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i)))−ln⁡p(Y∣θ(i))=∑zp(Z∣Y,θ(i))ln⁡(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))p(Y∣θ(i))).

于是,令B(θ,θ(i))=defH(θ(i)+∑zp(Z∣Y,θ(i))ln(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))p(Y∣θ(i))),B(θ,θ(i))=defH(θ(i)+∑zp(Z∣Y,θ(i))ln⁡(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))p(Y∣θ(i))),

则有 H(θ)⩾B(θ,θ(i))H(θ)⩾B(θ,θ(i)),并且有 H(θ(i))=B(θ(i),θ(i))H(θ(i))=B(θ(i),θ(i))。因此,实际上我们已经构造了 H(θ)H(θ) 的一个下界 B(θ,θ(i))B(θ,θ(i)),即完成了 E 步。接下来,我们要想办法提升这个下界(M 步)。

优化目标

考虑到提升 B(θ,θ(i))B(θ,θ(i)) 就能提升 H(θ)H(θ)。原始的优化问题转化为:θ(i+1)=====argmaxθB(θ,θ(i))argmaxθ[H(θ(i))+∑zp(Z∣Y,θ(i))ln(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))p(Y∣θ(i)))]argmaxθ[∑zp(Z∣Y,θ(i))ln(p(Z∣θ)p(Y∣Z,θ))+H(θ(i))−∑zp(Z∣Y,θ(i))ln(p(Z∣Y,θ(i))p(Y∣θ(i)))]argmaxθ[∑zp(Z∣Y,θ(i))ln(p(Z∣θ)p(Y∣Z,θ))]// 省略相对 θ 的常数项argmaxθ[∑zp(Z∣Y,θ(i))lnp(Y,Z∣θ)]θ(i+1)=argmaxθB(θ,θ(i))=argmaxθ[H(θ(i))+∑zp(Z∣Y,θ(i))ln⁡(p(Z∣θ)p(Y∣Z,θ)p(Z∣Y,θ(i))p(Y∣θ(i)))]=argmaxθ[∑zp(Z∣Y,θ(i))ln⁡(p(Z∣θ)p(Y∣Z,θ))+H(θ(i))−∑zp(Z∣Y,θ(i))ln⁡(p(Z∣Y,θ(i))p(Y∣θ(i)))]=argmaxθ[∑zp(Z∣Y,θ(i))ln⁡(p(Z∣θ)p(Y∣Z,θ))]// 省略相对 θ 的常数项=argmaxθ[∑zp(Z∣Y,θ(i))ln⁡p(Y,Z∣θ)]

于是,令Q(θ,θ(i))=def∑zp(Z∣Y,θ(i))lnp(Y,Z∣θ).Q(θ,θ(i))=def∑zp(Z∣Y,θ(i))ln⁡p(Y,Z∣θ).

这即是在 M 步我们需要优化的目标。

EM 迭代算法

  • 输入:观测变量数据 YY,隐变量数据 ZZ,给定参数 θθ 时观测变量和隐变量的联合分布 p(Y,Z∣θ)p(Y,Z∣θ),给定观测变量和参数 θθ 时隐变量 ZZ 的条件分布 p(Z∣Y,θ)p(Z∣Y,θ)。
  • 输出:模型参数 θθ。
  • 设置参数初值 θ(0)θ(0)。
  • 开始迭代,直至 θθ 收敛:
    • E 步,计算Q(θ,θ(i))=def∑zp(Z∣Y,θ(i))lnp(Y,Z∣θ).Q(θ,θ(i))=def∑zp(Z∣Y,θ(i))ln⁡p(Y,Z∣θ).
    • M 步,解最优化问题θ(i+1)=argmaxθQ(θ,θ(i)).θ(i+1)=argmaxθQ(θ,θ(i)).

可见,在模型输入中包含联合分布和条件分布。这印证了我们之前说的:「EM 算法同极大似然法一样,需要预先假定分布」。

EM 算法的收敛性

考虑 EM 算法在迭代过程中产生的参数估计序列 θ(i)θ(i),以及据此计算的对数似然函数序列 H(θ(i))=lnp(Y∣θ(i)=B(θ(i),θ(i))H(θ(i))=ln⁡p(Y∣θ(i)=B(θ(i),θ(i))。

考虑到在 EM 迭代算法中有:B(θ(i+1),θ(i+1))⩾B(θ(i+1),θ(i))⩾B(θ(i),θ(i)).B(θ(i+1),θ(i+1))⩾B(θ(i+1),θ(i))⩾B(θ(i),θ(i)).

因此,上述对数似然函数序列是单调递增的。又考虑到似然函数是概率,故而上述对数似然函数序列有上界。根据单调有界原理,上述对数似然函数序列必然收敛。不过,这并不意味着参数的估计序列 θ(i)θ(i) 必然收敛。因此,在实际使用过程中,EM 算法的终止条件除了 ∥∥θ(i+1)−θ(i)∥∥<ε1‖θ(i+1)−θ(i)‖<ε1 之外,还应包含 ∥∥Q(θ(i+1),θ(i))−Q(θ(i),θ(i))∥∥<ε2‖Q(θ(i+1),θ(i))−Q(θ(i),θ(i))‖<ε2。

此外,由于 EM 算法对初值敏感。故而,尽管 EM 算法能保证收敛,但不保证一定能收敛到全局最优解。因此,使用 EM 算法时,往往需要设置多个初值,估计多套模型参数。最后再将模型参数代入模型,比较模型的效果,选出其中最好的作为最终结果。

原文链接:https://liam.page/2018/11/08/Expectation-Maximization-Algorithm/

  • 微信(WeChat)
  • 多少不重要
  • weinxin
  • 支付宝(Alipay)
  • 有鼓励就好
  • weinxin
voice 站点
Combantrin澳大利亚驱虫巧克力安全高效美味成人宝宝打虫药
兰芝 多效洁面乳
进口洗发水一件代发,微商货源,韩国RYOE吕 洗发两支装 红吕
懒人床上笔记本电脑桌台式家用双人电脑桌床上书桌可移动跨床桌

发表评论

:?::razz::sad::evil::!::smile::oops::grin::eek::shock::???::cool::lol::mad::twisted::roll::wink::idea::arrow::neutral::cry::mrgreen: