EM算法

来自集智百科 - 复杂系统|人工智能|复杂科学|复杂网络|自组织
薄荷（讨论 | 贡献）2021年12月28日 (二) 19:59的版本 （→‎编者推荐）
跳到导航 跳到搜索


在统计学中，期望最大化算法 expectation–maximization algorithm（EM algorithm）是一种寻找统计模型中（局部）极大似然或者最大后验 maximum a posteriori（MAP）参数估计的迭代方法，其中的统计模型依赖于未观测到的潜在变量。EM迭代过程中交替执行期望（E）步和最大化（M）步；前者使用当前参数估计值建立对数似然函数的期望函数，后者计算能够最大化E步中获得的期望对数似然函数的参数。这些参数估计值将用于确定下一个E步中潜在变量的分布。

Old Faithful火山喷发数据的 EM 聚类。随机初始模型（由于轴的不同尺度，看起来是两个非常平和宽的区域）可以拟合观测的数据。在第一次迭代中，模型发生了实质性的变化，但随后收敛到间歇泉的两个模态。可视化使用 ELKI.


历史

Arthur Dempster, Nan Laird和Donald Rubin[1]于1977年发表的一篇经典的论文中解释和命名了EM算法。他们指出该方法被早期作者“多次在特殊条件下提出”。Cedric Smith提出的估计等位基因频率的基因计数法便是其中之一。[2]基于与Per Martin-Löf和Anders Martin-Löf的合作，Rolf Sundberg在他的学位论文和若干论文[3][4][5]中详述了针对指数族的EM方法。[6][7][8][9][10]1977年Dempster–Laird–Rubin的论文推广了该方法并针对更为广泛的问题进行了收敛性分析，该论文确立EM方法作为统计分析的重要工具。[11][12] Dempster–Laird–Rubin算法的收敛分析存在缺陷，C. F. Jeff Wu在1983年发表了一项修正的收敛性分析。[13]Dempster–Laird–Rubin称，Wu的工作建立了EM方法在指数族之外的收敛。[13]


介绍

EM 算法用于在无法直接求解方程的情况下查找统计模型的（局部）最大似然参数。通常，除了未知参数和已知数据观察之外，这些模型还涉及潜在变量。也就是说，要么数据中存在缺失值，要么可以通过假设存在更多未观察到的数据点来更简单地制定模型。例如，通过假设每个观察到的数据点都有一个对应的未观察到的数据点或潜在变量，指定每个数据点所属的混合成分，可以更简单地描述混合模型。


寻找最大似然解通常需要对所有未知值、参数和潜在变量求似然函数的导数，并同时求解所得方程。在具有潜在变量的统计模型中，这通常是不可能的。相反，结果通常是一组互锁方程，其中参数的解需要潜在变量的值，反之亦然，但将一组方程代入另一组会产生无法求解的方程。


EM 算法从观察开始到找到方法可以数值求解这两组方程。可以简单地为两组未知数中的一组选择任意值，使用它们来估计第二组，然后使用这些新值来找到对第一组更好的估计，然后在两者之间保持交替，直到结果值都收敛到固定点。[13] 这并不明显，但可以证明在这种情况下它确实有效，并且可能性的导数在该点（任意接近）为零，这反过来意味着该点是最大值或一个鞍点。通常，可能会出现多个最大值，但不能保证会找到全局最大值。一些可能性也有奇点，即无意义的最大值。例如，EM 在混合模型中可能找到的解决方案之一包含把其中一个成分设置为具有零方差，并且相同成分的平均参数和一个数据点等价。


描述

给定生成一组 [math]\displaystyle{ \mathbf{X} }[/math] 观察数据，一组未观察到的潜在数据或缺失值 [math]\displaystyle{ \mathbf{Z} }[/math]的统计模型， 和未知参数向量 [math]\displaystyle{ \boldsymbol\theta }[/math]，以及似然函数 [math]\displaystyle{ L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z}) = p(\mathbf{X}, \mathbf{Z}\mid\boldsymbol\theta) }[/math]，未知参数的最大似然估计 maximum likelihood estimate（MLE） 通过最大化观察数据的边际可能性来决定。


[math]\displaystyle{ L(\boldsymbol\theta; \mathbf{X}) = p(\mathbf{X}\mid\boldsymbol\theta) = \int p(\mathbf{X},\mathbf{Z} \mid \boldsymbol\theta) \, d\mathbf{Z} }[/math]

然而，这个量通常是难以处理的，因为[math]\displaystyle{ \mathbf{Z} }[/math]是不可观察的，并且[math]\displaystyle{ \mathbf{Z} }[/math] 的分布在达到[math]\displaystyle{ \boldsymbol\theta }[/math]之前是未知的。


EM算法试图通过迭代应用这两个步骤来找到边际可能性的最大似然估计：

期望步（E步）: 定义[math]\displaystyle{ Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) }[/math]作为 [math]\displaystyle{ \boldsymbol\theta }[/math] 的对数似然函数的期望值 ，关于 [math]\displaystyle{ \mathbf{Z} }[/math] 的当前条件分布给定 [math]\displaystyle{ \mathbf{X} }[/math] 和参数[math]\displaystyle{ \boldsymbol\theta^{(t)} }[/math]:


[math]\displaystyle{ Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) = \operatorname{E}_{\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}}\left[ \log L (\boldsymbol\theta; \mathbf{X},\mathbf{Z}) \right] \, }[/math]
最大化步（M步）: 找到使数量最大化的参数：
[math]\displaystyle{ \boldsymbol\theta^{(t+1)} = \underset{\boldsymbol\theta}{\operatorname{arg\,max}} \ Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) \, }[/math]
应用 EM 的典型模型使用[math]\displaystyle{ \mathbf{Z} }[/math]作为潜在变量，表明属于一组组之一：
1. 观察到的数据点[math]\displaystyle{ \mathbf{X} }[/math]可能是离散的（取有限或可数无限集合中的值）或连续（取不可数无限集合中的值）。 与每个数据点相关联的可能是观察向量。
2. 缺失值（又名潜在变量）[math]\displaystyle{ \mathbf{Z} }[/math]是离散的，从固定数量的值中提取，每个观察单位有一个潜在变量。
3. 参数是连续的，分为两种：与所有数据点相关的参数，以及与潜在变量的特定值相关的参数（即与对应潜在变量具有该值的所有数据点相关联）。


然而，EM也可以应用到其他模型中。

原因如下。如果已知参数 [math]\displaystyle{ \boldsymbol\theta }[/math]的值，通常可以通过最大化 [math]\displaystyle{ \mathbf{Z} }[/math] 的所有可能值的对数似然找到潜变量[math]\displaystyle{ \mathbf{Z} }[/math]的值，或者简单地通过迭代 [math]\displaystyle{ \mathbf{Z} }[/math]或通过一种算法，例如用于隐马尔可夫模型的 Baum-Welch 算法。相反，如果我们知道潜在变量[math]\displaystyle{ \mathbf{Z} }[/math]的值，我们可以非常容易找到参数 [math]\displaystyle{ \boldsymbol\theta }[/math] 的估计，通常只需根据相关潜在变量的值对观察到的数据点进行分组，并对每组中点的值或值的某些函数求平均值。这表明了一种在[math]\displaystyle{ \boldsymbol\theta }[/math]和 [math]\displaystyle{ \mathbf{Z} }[/math] 都是未知情况下的迭代算法：


  1. 首先，将参数[math]\displaystyle{ \boldsymbol\theta }[/math] 初始化为一些随机值。
  2. 计算[math]\displaystyle{ \mathbf{Z} }[/math] 的每个可能值的概率，给定[math]\displaystyle{ \boldsymbol\theta }[/math]。
  3. 然后，使用刚刚计算的 [math]\displaystyle{ \mathbf{Z} }[/math]来计算参数[math]\displaystyle{ \boldsymbol\theta }[/math]的更好估计。
  4. 迭代步骤 2 和 3 直到收敛。


刚刚描述的算法单调地接近成本函数的局部最小值。


特性

说到期望 (E) 步骤有点用词不当。 在第一步中计算的是函数Q的固定的、与数据相关的参数。一旦Q的参数已知，就可以完全确定并在 EM 算法的第二步 (M) 中最大化。


尽管 EM 迭代确实增加了观察到的数据（即边际）似然函数，但不能保证序列收敛到最大似然估计量。 对于多峰分布，这意味着 EM 算法可能会收敛到观察到的数据似然函数的局部最大值，这取决于起始值。 存在各种启发式或元启发式方法来逃避局部最大值，例如随机重启爬山（从几个不同的随机初始估计值θ(t)开始），或应用模拟退火方法。


当似然是指数族时，EM 尤其有用：E 步成为充分统计量的期望总和，而 M 步涉及最大化线性函数。 在这种情况下，通常可以使用 Sundberg 公式（由 Rolf Sundberg 发布，使用 Per Martin-Löf 和 Anders Martin-Löf 的未发表结果）为每个步骤导出封闭形式的表达式更新。


在 Dempster、Laird 和 Rubin 的原始论文中，EM 方法被修改为计算贝叶斯推理的最大后验估计。


存在其他方法来寻找最大似然估计，例如梯度下降、共轭梯度或高斯-牛顿算法的变体。 与 EM 不同，此类方法通常需要评估似然函数的一阶和/或二阶导数。


正确性证明

期望最大化可以改善[math]\displaystyle{ Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) }[/math]而不是直接改进[math]\displaystyle{ \log p(\mathbf{X}\mid\boldsymbol\theta) }[/math] 。 这里表明对前者的改进意味着对后者的改进。[14]


对于任何具有非零概率 [math]\displaystyle{ \mathbf{Z} }[/math] with non-zero probability [math]\displaystyle{ p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta) }[/math]，我们可以写

[math]\displaystyle{ \log p(\mathbf{X}\mid\boldsymbol\theta) = \log p(\mathbf{X},\mathbf{Z}\mid\boldsymbol\theta) - \log p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta). }[/math]


我们在当前参数估计[math]\displaystyle{ \theta^{(t)} }[/math]下对未知数据的可能值取期望值[math]\displaystyle{ \mathbf{Z} }[/math]两边乘以[math]\displaystyle{ p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}) }[/math] 并在 [math]\displaystyle{ \mathbf{Z} }[/math]上求和（或积分）。左边是一个常数的期望，所以我们得到：

[math]\displaystyle{ \begin{align} \log p(\mathbf{X}\mid\boldsymbol\theta) & = \sum_{\mathbf{Z}} p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}) \log p(\mathbf{X},\mathbf{Z}\mid\boldsymbol\theta) - \sum_{\mathbf{Z}} p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta^{(t)}) \log p(\mathbf{Z}\mid\mathbf{X},\boldsymbol\theta) \\ & = Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}), \end{align} }[/math]


其中 [math]\displaystyle{ H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) }[/math] 由它正在替换的否定和定义。最后一个方程适用于[math]\displaystyle{ \boldsymbol\theta }[/math] 的每个值，包括 [math]\displaystyle{ \boldsymbol\theta = \boldsymbol\theta^{(t)} }[/math],

[math]\displaystyle{ \log p(\mathbf{X}\mid\boldsymbol\theta^{(t)}) = Q(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}), }[/math]


并从前一个方程中减去最后一个方程给出

[math]\displaystyle{ \log p(\mathbf{X}\mid\boldsymbol\theta) - \log p(\mathbf{X}\mid\boldsymbol\theta^{(t)}) = Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) - Q(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}) + H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) - H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}), }[/math]


然而，吉布斯不等式告诉我们 [math]\displaystyle{ H(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) \ge H(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}) }[/math]，所以我们可以得出结论

[math]\displaystyle{ \log p(\mathbf{X}\mid\boldsymbol\theta) - \log p(\mathbf{X}\mid\boldsymbol\theta^{(t)}) \ge Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) - Q(\boldsymbol\theta^{(t)}\mid\boldsymbol\theta^{(t)}). }[/math]


换句话说，选择[math]\displaystyle{ \boldsymbol\theta }[/math]来改进[math]\displaystyle{ Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) }[/math] 导致 [math]\displaystyle{ \log p(\mathbf{X}\mid\boldsymbol\theta) }[/math] 至少有同样的改进。


作为最大化-最大化过程

EM算法可以看作是两个交替的最大化步骤，即作为坐标下降法的一个例子。[15][16]考虑函数：

[math]\displaystyle{ F(q,\theta) := \operatorname{E}_q [ \log L (\theta ; x,Z) ] + H(q), }[/math]


q是未观测数据 z上的任意概率分布，H(q)是分布q的熵。 这个函数可以写成

[math]\displaystyle{ F(q,\theta) = -D_{\mathrm{KL}}\big(q \parallel p_{Z\mid X}(\cdot\mid x;\theta ) \big) + \log L(\theta;x), }[/math]


其中[math]\displaystyle{ p_{Z\mid X}(\cdot\mid x;\theta ) }[/math]是在给定观察数据[math]\displaystyle{ x }[/math]前提下未观察到数据的条件分布，[math]\displaystyle{ D_{KL} }[/math]是 Kullback–Leibler 散度。那么EM算法的步骤可以看成：


期望步 Expectation step：选择[math]\displaystyle{ q }[/math]最大化[math]\displaystyle{ F }[/math]：

[math]\displaystyle{ q^{(t)} = \operatorname{arg\,max}_q \ F(q,\theta^{(t)}) }[/math]


最大化步 Maximization step: 选择[math]\displaystyle{ \theta }[/math]最大化[math]\displaystyle{ F }[/math]:


应用

EM 经常用于机器学习和计算机视觉中的数据聚类。 在自然语言处理中，该算法的两个突出实例是用于隐马尔可夫模型的 Baum-Welch 算法和用于无监督概率上下文无关文法归纳的内-外算法。EM 经常用于混合模型的参数估计，[17][18]特别是在数量遗传学中。[19]


在心理测量学中，EM对于估计项目响应理论模型的项目参数和潜在能力几乎是必不可少的。


由于能够处理丢失的数据和观察不明变量，EM 正成为一个有用的对投资组合进行定价和管理风险的工具。


EM 算法(及其更快变体的有序子集期望最大化)也广泛应用于医学图像重建，特别是在正电子发射断层扫描、单光子发射计算机断层扫描和X射线计算机断层扫描中。下面是 EM 的其他更快的变体。


在结构工程中，利用期望最大化(STRIDE)[20]算法进行结构识别是一种仅输出的方法，用于利用传感器数据识别结构系统的自然振动特性(见运行模态分析)。


滤波和平滑EM算法

卡尔曼滤波器通常用于在线状态估计，最小方差平滑器可用于离线或批状态估计。然而，这些最小方差解需要状态空间模型参数的估计。EM 算法可用于求解联合状态和参数估计问题。


滤波和平滑 EM 算法是通过重复这两个步骤产生的:


E步

操作一个 Kalman 滤波器或一个最小方差平滑设计与当前的参数估计，以获得更新的状态估计。


M步

使用最大似然计算中的滤波或平滑状态估计来获得更新的参数估计。


假设卡尔曼滤波器或最小方差平滑器对具有附加白噪声的单输入单输出系统进行测量。通过极大似然估计，可以得到更新的测量噪声方差估计

[math]\displaystyle{ \widehat{\sigma}^2_v = \frac{1}{N} \sum_{k=1}^N {(z_k-\widehat{x}_k)}^2, }[/math]


其中[math]\displaystyle{ \widehat{x}_k }[/math]是由滤波器或平滑器从 N 个标量测量[math]\displaystyle{ z_k }[/math]。 上述更新也可以应用于更新泊松测量噪声强度。 类似地，对于一阶自回归过程，更新的过程噪声方差估计可以计算为</nowiki>

[math]\displaystyle{ \widehat{\sigma}^2_w = \frac{1}{N} \sum_{k=1}^N {(\widehat{x}_{k+1}-\widehat{F}\widehat{{x}}_k)}^2, }[/math]


其中 [math]\displaystyle{ \widehat{x}_k }[/math] 和 [math]\displaystyle{ \widehat{x}_{k+1} }[/math]是由过滤器或平滑器计算的标量状态估计。 更新后的模型系数估计是通过</nowiki>

[math]\displaystyle{ \widehat{F} = \frac{\sum_{k=1}^N (\widehat{x}_{k+1}-\widehat{F} \widehat{x}_k)}{\sum_{k=1}^N \widehat{x}_k^2}. }[/math]


研究了上述参数估计的收敛性。[21][22][23][24]


变体

针对有时EM 算法收敛速度慢的问题，一些改进方法被提出，如共轭梯度法和修正牛顿法(Newton-Raphson)。[25] 此外，EM 还可以与约束估计方法一起使用。


参数扩展期望最大化 Parameter-expanded expectation maximization (PX-EM)算法通过“利用输入完整数据中获得的额外信息，通过‘协方差调整’来校正 m 步的分析” ，从而提高了计算速度。[26]


期望条件最大化 Expectation conditional maximization (ECM)用一系列条件最大化(CM)步骤代替每个 m 步骤，其中每个参数θi单独最大化，条件是其他参数保持不变。[27] 本身也可以扩展为期望条件最大化 Expectation conditional maximization either(ECME)算法。[28]


也可以考虑将 EM 算法作为 MM (majorize/minize 或 Minorize/Maximize，取决于上下文)算法的子类，[29]因此可以使用在更一般情况下开发的任何机制。


α-EM算法

EM 算法中使用的 Q 函数基于对数似然。 因此，它被视为log-EM算法。 对数似然的使用可以推广到 α-对数似然比的使用。 然后，通过使用α-log似然比和α-散度的Q函数，可以将观测数据的α-log似然比精确表示为等式。 获得这个 Q 函数是一个广义的 E 步骤。 它的最大化是一个广义的 M 步。 这对称为 α-EM 算法，[30]它包含 log-EM 算法作为其子类。 因此，Yasuo Matsuyama 的 α-EM 算法是 log-EM 算法的精确推广。 不需要计算梯度或 Hessian 矩阵。 通过选择合适的 α，α-EM 显示出比 log-EM 算法更快的收敛速度。 α-EM 算法导致了隐马尔可夫模型估计算法 α-HMM 的更快版本。[31]


与变分贝叶斯方法的关系

EM 是一个部分非贝叶斯，最大似然方法。它的最终结果给出了一个关于潜在变量的概率分布估计(在贝叶斯风格)以及 θ 的点估计(无论是最大似然估计还是后验模式)。一个完整的贝叶斯版本即给出一个关于θ 和潜在变量的概率分布可能是想要的。贝叶斯推理方法简单地将 θ 作为另一个潜变量来处理。在这个范例中，E 和 M 步骤之间的区别就消失了。如果使用上述因子化 Q 近似(变分贝叶斯) ，求解可以迭代每个潜变量(现在包括 θ) 并每次优化一个。现在，每次迭代需要 k 个步骤，其中 k 是潜变量的数量。对于图形模型，这是很容易做到的，因为每个变量的新 Q 只依赖于它的马尔可夫包层，所以局部信息传递可以用于有效的推理。


几何解释

在信息几何中，E步和M步被解释为双重仿射连接下的投影，称为e-connection和m-connection； Kullback-Leibler 背离也可以用这些术语来理解。


例子

让[math]\displaystyle{ \mathbf{x} = (\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_n) }[/math] 成为一个样本 [math]\displaystyle{ n }[/math]来自维度的两个多元正态分布的混合的独立观察[math]\displaystyle{ d }[/math],然后让[math]\displaystyle{ \mathbf{z} = (z_1,z_2,\ldots,z_n) }[/math]是确定观察来源的成分的潜在变量。[16]

[math]\displaystyle{ X_i \mid(Z_i = 1) \sim \mathcal{N}_d(\boldsymbol{\mu}_1,\Sigma_1) }[/math] and [math]\displaystyle{ X_i \mid(Z_i = 2) \sim \mathcal{N}_d(\boldsymbol{\mu}_2,\Sigma_2), }[/math]

其中

[math]\displaystyle{ \operatorname{P} (Z_i = 1 ) = \tau_1 \, }[/math] and [math]\displaystyle{ \operatorname{P} (Z_i=2) = \tau_2 = 1-\tau_1. }[/math]

目的是估计代表高斯分布之间混合值的未知参数以及每个的均值和协方差：

[math]\displaystyle{ \theta = \big( \boldsymbol{\tau},\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\Sigma_1,\Sigma_2 \big), }[/math]


其中不完全数据似然函数是

[math]\displaystyle{ L(\theta;\mathbf{x}) = \prod_{i=1}^n \sum_{j=1}^2 \tau_j \ f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j), }[/math]

并且完全数据似然函数是

[math]\displaystyle{ L(\theta;\mathbf{x},\mathbf{z}) = p(\mathbf{x},\mathbf{z} \mid \theta) = \prod_{i=1}^n \prod_{j=1}^2 \ [f(\mathbf{x}_i;\boldsymbol{\mu}_j,\Sigma_j) \tau_j] ^{\mathbb{I}(z_i=j)}, }[/math]

或者

[math]\displaystyle{ L(\theta;\mathbf{x},\mathbf{z}) = \exp \left\{ \sum_{i=1}^n \sum_{j=1}^2 \mathbb{I}(z_i=j) \big[ \log \tau_j -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{d}{2} \log(2\pi) \big] \right\}, }[/math]


其中 [math]\displaystyle{ \mathbb{I} }[/math]是一个指示函数，并且[math]\displaystyle{ f }[/math]是多元正态的概率密度函数。


在最后一个等式中，对于每个i，一个指标 [math]\displaystyle{ \mathbb{I}(z_i=j) }[/math]等于零，一个指标等于一。因此，内和减少为一项。


E 步骤

鉴于我们当前对参数θ(t)估计， Zi的条件分布由贝叶斯定理确定为由 τ加权的正态密度的比例高度：

[math]\displaystyle{ T_{j,i}^{(t)} := \operatorname{P}(Z_i=j \mid X_i=\mathbf{x}_i ;\theta^{(t)}) = \frac{\tau_j^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_j^{(t)},\Sigma_j^{(t)})}{\tau_1^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_1^{(t)},\Sigma_1^{(t)}) + \tau_2^{(t)} \ f(\mathbf{x}_i;\boldsymbol{\mu}_2^{(t)},\Sigma_2^{(t)})}. }[/math]

这些被称为“成员概率”，通常被认为是 E 步骤的输出（尽管这不是下面的 Q 函数）。


此 E 步骤对应于为 Q 设置此功能：

[math]\displaystyle{ \begin{align}Q(\theta\mid\theta^{(t)}) &= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X},\mathbf{\theta}^{(t)}} [\log L(\theta;\mathbf{x},\mathbf{Z}) ] \\ &= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X},\mathbf{\theta}^{(t)}} [\log \prod_{i=1}^{n}L(\theta;\mathbf{x}_i,Z_i) ] \\ &= \operatorname{E}_{\mathbf{Z}\mid\mathbf{X},\mathbf{\theta}^{(t)}} [\sum_{i=1}^n \log L(\theta;\mathbf{x}_i,Z_i) ] \\ &= \sum_{i=1}^n\operatorname{E}_{Z_i\mid\mathbf{X};\mathbf{\theta}^{(t)}} [\log L(\theta;\mathbf{x}_i,Z_i) ] \\ &= \sum_{i=1}^n \sum_{j=1}^2 P(Z_i =j \mid X_i = \mathbf{x}_i; \theta^{(t)}) \log L(\theta_j;\mathbf{x}_i,j) \\ &= \sum_{i=1}^n \sum_{j=1}^2 T_{j,i}^{(t)} \big[ \log \tau_j -\tfrac{1}{2} \log |\Sigma_j| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_j)^\top\Sigma_j^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_j) -\tfrac{d}{2} \log(2\pi) \big]. \end{align} }[/math]

[math]\displaystyle{ \log L(\theta;\mathbf{x}_i,Z_i) }[/math]的期望求和的内部是[math]\displaystyle{ P(Z_i \mid X_i = \mathbf{x}_i; \theta^{(t)}) }[/math]的概率密度函数这可能因人而异[math]\displaystyle{ \mathbf{x}_i }[/math] 的训练集。E 步骤中的所有内容在执行该步骤之前都是已知的，除了根据 E 步骤部分开头的方程计算的[math]\displaystyle{ T_{j,i} }[/math]。


这一完整的条件期望不需要一步计算，因为 τ和 μ/Σ 出现在单独的线性项中，因此可以独立最大化。


M 步骤

Q(θ | θ(t))在形式上是二次的，这意味着确定 θ 最大值相对简单。此外， τ, (μ1,Σ1)和(μ2,Σ2)都可以独立最大化，因为它们都出现在单独的线性项中。


首先，考虑τ，其具有约束 τ1 + τ2=1:

[math]\displaystyle{ \begin{align}\boldsymbol{\tau}^{(t+1)} &= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}}\ Q(\theta \mid \theta^{(t)} ) \\ &= \underset{\boldsymbol{\tau}} {\operatorname{arg\,max}} \ \left\{ \left[ \sum_{i=1}^n T_{1,i}^{(t)} \right] \log \tau_1 + \left[ \sum_{i=1}^n T_{2,i}^{(t)} \right] \log \tau_2 \right\}. \end{align} }[/math]


这与二项式分布的 MLE 具有相同的形式，因此

[math]\displaystyle{ \tau^{(t+1)}_j = \frac{\sum_{i=1}^n T_{j,i}^{(t)}}{\sum_{i=1}^n (T_{1,i}^{(t)} + T_{2,i}^{(t)} ) } = \frac{1}{n} \sum_{i=1}^n T_{j,i}^{(t)}. }[/math]


对于 (μ1,Σ1)的下一个估计：

[math]\displaystyle{ \begin{align}(\boldsymbol{\mu}_1^{(t+1)},\Sigma_1^{(t+1)}) &= \underset{\boldsymbol{\mu}_1,\Sigma_1} {argmax} Q(\theta \mid \theta^{(t)} ) \\ &= \underset{\boldsymbol{\mu}_1,\Sigma_1} {argmax} \sum_{i=1}^n T_{1,i}^{(t)} \left\{ -\tfrac{1}{2} \log |\Sigma_1| -\tfrac{1}{2}(\mathbf{x}_i-\boldsymbol{\mu}_1)^\top\Sigma_1^{-1} (\mathbf{x}_i-\boldsymbol{\mu}_1) \right\} \end{align}. }[/math]


这与正态分布的加权 MLE 具有相同的形式，因此

[math]\displaystyle{ \boldsymbol{\mu}_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} \mathbf{x}_i}{\sum_{i=1}^n T_{1,i}^{(t)}} }[/math] and [math]\displaystyle{ \Sigma_1^{(t+1)} = \frac{\sum_{i=1}^n T_{1,i}^{(t)} (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)}) (\mathbf{x}_i - \boldsymbol{\mu}_1^{(t+1)})^\top }{\sum_{i=1}^n T_{1,i}^{(t)}} }[/math]


并且，通过对称性，

[math]\displaystyle{ \boldsymbol{\mu}_2^{(t+1)} = \frac{\sum_{i=1}^n T_{2,i}^{(t)} \mathbf{x}_i}{\sum_{i=1}^n T_{2,i}^{(t)}} }[/math] and [math]\displaystyle{ \Sigma_2^{(t+1)} = \frac{\sum_{i=1}^n T_{2,i}^{(t)} (\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)}) (\mathbf{x}_i - \boldsymbol{\mu}_2^{(t+1)})^\top }{\sum_{i=1}^n T_{2,i}^{(t)}}。 }[/math]


终止

如果[math]\displaystyle{ E_{Z\mid\theta^{(t)},\mathbf{x}}[\log L(\theta^{(t)};\mathbf{x},\mathbf{Z})] \leq E_{Z\mid\theta^{(t-1)},\mathbf{x}}[\log L(\theta^{(t-1)};\mathbf{x},\mathbf{Z})] + \varepsilon }[/math] for [math]\displaystyle{ \varepsilon }[/math] 低于某个预设阈值终止迭代过程。

一般化

上面说明的算法可以推广到两个以上多元正态分布的混合。


截断和删减回归

EM 算法已在存在解释某些量变化的基础线性回归模型的情况下实施，但实际观察到的值是模型中表示的那些值的删失或截断版本。 此模型的特殊情况包括来自一个正态分布的删失或截断观察。


参考文献

  1. ↑ Dempster, A.P.; Laird, N.M.; Rubin, D.B. (1977). "Maximum Likelihood from Incomplete Data via the EM Algorithm". Journal of the Royal Statistical Society, Series B. 39 (1): 1–38. JSTOR 2984875. MR 0501537.
  2. ↑ Ceppelini, R.M. (1955). "The estimation of gene frequencies in a random-mating population". Ann. Hum. Genet. 20 (2): 97–115. doi:10.1111/j.1469-1809.1955.tb01360.x. PMID 13268982.
  3. ↑ Sundberg, Rolf (1974). "Maximum likelihood theory for incomplete data from an exponential family". Scandinavian Journal of Statistics. 1 (2): 49–58. JSTOR 4615553. MR 0381110.
  4. ↑ Rolf Sundberg. 1971. Maximum likelihood theory and applications for distributions generated when observing a function of an exponential family variable. Dissertation, Institute for Mathematical Statistics, Stockholm University.
  5. ↑ Sundberg, Rolf (1976). "An iterative method for solution of the likelihood equations for incomplete data from exponential families". Communications in Statistics – Simulation and Computation. 5 (1): 55–64. doi:10.1080/03610917608812007. MR 0443190.
  6. ↑ See the acknowledgement by Dempster, Laird and Rubin on pages 3, 5 and 11.
  7. ↑ G. Kulldorff. 1961. Contributions to the theory of estimation from grouped and partially grouped samples. Almqvist & Wiksell.
  8. ↑ Anders Martin-Löf. 1963. "Utvärdering av livslängder i subnanosekundsområdet" ("Evaluation of sub-nanosecond lifetimes"). ("Sundberg formula")
  9. ↑ Per Martin-Löf. 1966. Statistics from the point of view of statistical mechanics. Lecture notes, Mathematical Institute, Aarhus University. ("Sundberg formula" credited to Anders Martin-Löf).
  10. ↑ Per Martin-Löf. 1970. Statistika Modeller (Statistical Models): Anteckningar från seminarier läsåret 1969–1970 (Notes from seminars in the academic year 1969-1970), with the assistance of Rolf Sundberg. Stockholm University. ("Sundberg formula")
  11. ↑ Martin-Löf, P. The notion of redundancy and its use as a quantitative measure of the deviation between a statistical hypothesis and a set of observational data. With a discussion by F. Abildgård, A. P. Dempster, D. Basu, D. R. Cox, A. W. F. Edwards, D. A. Sprott, G. A. Barnard, O. Barndorff-Nielsen, J. D. Kalbfleisch and G. Rasch and a reply by the author. Proceedings of Conference on Foundational Questions in Statistical Inference (Aarhus, 1973), pp. 1–42. Memoirs, No. 1, Dept. Theoret. Statist., Inst. Math., Univ. Aarhus, Aarhus, 1974.
  12. ↑ Martin-Löf, Per (1974). "The notion of redundancy and its use as a quantitative measure of the discrepancy between a statistical hypothesis and a set of observational data". Scand. J. Statist. 1 (1): 3–18.
  13. ↑ 13.0 13.1 13.2 Wu, C. F. Jeff (Mar 1983). "On the Convergence Properties of the EM Algorithm". Annals of Statistics. 11 (1): 95–103. doi:10.1214/aos/1176346060. JSTOR 2240463. MR 0684867.
  14. ↑ Little, Roderick J.A.; Rubin, Donald B. (1987). Statistical Analysis with Missing Data. Wiley Series in Probability and Mathematical Statistics. New York: John Wiley & Sons. pp. 134–136. ISBN 978-0-471-80254-9. https://archive.org/details/statisticalanaly00litt. 
  15. ↑ Neal, Radford; Hinton, Geoffrey (1999). Michael I. Jordan. ed. A view of the EM algorithm that justifies incremental, sparse, and other variants. Cambridge, MA: MIT Press. pp. 355–368. ISBN 978-0-262-60032-3. http://ftp.cs.toronto.edu/pub/radford/emk.pdf. 
  16. ↑ 16.0 16.1 Hastie, Trevor; Tibshirani, Robert; Friedman, Jerome (2001). "8.5 The EM algorithm". The Elements of Statistical Learning. New York: Springer. pp. 236–243. ISBN 978-0-387-95284-0. https://archive.org/details/elementsstatisti00thas_842. 
  17. ↑ Lindstrom, Mary J; Bates, Douglas M (1988). "Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data". Journal of the American Statistical Association. 83 (404): 1014. doi:10.1080/01621459.1988.10478693.
  18. ↑ Van Dyk, David A (2000). "Fitting Mixed-Effects Models Using Efficient EM-Type Algorithms". Journal of Computational and Graphical Statistics. 9 (1): 78–98. doi:10.2307/1390614. JSTOR 1390614.
  19. ↑ Diffey, S. M; Smith, A. B; Welsh, A. H; Cullis, B. R (2017). "A new REML (parameter expanded) EM algorithm for linear mixed models". Australian & New Zealand Journal of Statistics. 59 (4): 433. doi:10.1111/anzs.12208.
  20. ↑ Matarazzo, T. J., and Pakzad, S. N. (2016). “STRIDE for Structural Identification using Expectation Maximization: Iterative Output-Only Method for Modal Identification.” Journal of Engineering Mechanics.http://ascelibrary.org/doi/abs/10.1061/(ASCE)EM.1943-7889.0000951
  21. ↑ Einicke, G. A.; Malos, J. T.; Reid, D. C.; Hainsworth, D. W. (January 2009). "Riccati Equation and EM Algorithm Convergence for Inertial Navigation Alignment". IEEE Trans. Signal Process. 57 (1): 370–375. Bibcode:2009ITSP...57..370E. doi:10.1109/TSP.2008.2007090.
  22. ↑ Einicke, G. A.; Falco, G.; Malos, J. T. (May 2010). "EM Algorithm State Matrix Estimation for Navigation". IEEE Signal Processing Letters. 17 (5): 437–440. Bibcode:2010ISPL...17..437E. doi:10.1109/LSP.2010.2043151.
  23. ↑ Einicke, G. A.; Falco, G.; Dunn, M. T.; Reid, D. C. (May 2012). "Iterative Smoother-Based Variance Estimation". IEEE Signal Processing Letters. 19 (5): 275–278. Bibcode:2012ISPL...19..275E. doi:10.1109/LSP.2012.2190278.
  24. ↑ Einicke, G. A. (Sep 2015). "Iterative Filtering and Smoothing of Measurements Possessing Poisson Noise". IEEE Transactions on Aerospace and Electronic Systems. 51 (3): 2205–2011. Bibcode:2015ITAES..51.2205E. doi:10.1109/TAES.2015.140843.
  25. ↑ Jamshidian, Mortaza; Jennrich, Robert I. (1997). "Acceleration of the EM Algorithm by using Quasi-Newton Methods". Journal of the Royal Statistical Society, Series B. 59 (2): 569–587. doi:10.1111/1467-9868.00083. MR 1452026.
  26. ↑ Liu, C (1998). "Parameter expansion to accelerate EM: The PX-EM algorithm". Biometrika. 85 (4): 755–770. CiteSeerX 10.1.1.134.9617. doi:10.1093/biomet/85.4.755.
  27. ↑ Meng, Xiao-Li; Rubin, Donald B. (1993). "Maximum likelihood estimation via the ECM algorithm: A general framework". Biometrika. 80 (2): 267–278. doi:10.1093/biomet/80.2.267. MR 1243503.
  28. ↑ Liu, Chuanhai; Rubin, Donald B (1994). "The ECME Algorithm: A Simple Extension of EM and ECM with Faster Monotone Convergence". Biometrika. 81 (4): 633. doi:10.1093/biomet/81.4.633. JSTOR 2337067.
  29. ↑ Hunter DR and Lange K (2004), A Tutorial on MM Algorithms, The American Statistician, 58: 30–37
  30. ↑ Matsuyama, Yasuo (2003). "The α-EM algorithm: Surrogate likelihood maximization using α-logarithmic information measures". IEEE Transactions on Information Theory. 49 (3): 692–706. doi:10.1109/TIT.2002.808105.
  31. ↑ Matsuyama, Yasuo (2011). "Hidden Markov model estimation based on alpha-EM algorithm: Discrete and continuous alpha-HMMs". International Joint Conference on Neural Networks: 808–816.


进一步阅读

  • Hogg, Robert; McKean, Joseph; Craig, Allen (2005). Introduction to Mathematical Statistics. Upper Saddle River, NJ: Pearson Prentice Hall. pp. 359–364. 
  • Dellaert, Frank (2002). "The Expectation Maximization Algorithm". CiteSeerX 10.1.1.9.9735. {{cite journal}}: Cite journal requires |journal= (help) gives an easier explanation of EM algorithm as to lowerbound maximization.
  • Bishop, Christopher M. (2006). Pattern Recognition and Machine Learning. Springer. ISBN 978-0-387-31073-2. 
  • Gupta, M. R.; Chen, Y. (2010). "Theory and Use of the EM Algorithm". Foundations and Trends in Signal Processing. 4 (3): 223–296. CiteSeerX 10.1.1.219.6830. doi:10.1561/2000000034. A well-written short book on EM, including detailed derivation of EM for GMMs, HMMs, and Dirichlet.
  • Bilmes, Jeff (1998). "A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models". CiteSeerX 10.1.1.28.613. {{cite journal}}: Cite journal requires |journal= (help) includes a simplified derivation of the EM equations for Gaussian Mixtures and Gaussian Mixture Hidden Markov Models.
  • McLachlan, Geoffrey J.; Krishnan, Thriyambakam (2008). The EM Algorithm and Extensions (2nd ed.). Hoboken: Wiley. ISBN 978-0-471-20170-0. 


编者推荐

相关文章

下一秒火焰如何燃烧？机器学习成功预测混沌

该文章是对Machine Learning’s ‘Amazing’ Ability to Predict Chaos的概述，介绍了机器学习预测混沌的理论基础及步骤。


论文速递集合


集智相关课程

课程：当复杂网络遇上机器学习

机器学习入门

该课程由莫烦、张江和尹相志三位老师共同教授，主要介绍了机器学习的相关理论及常用分析方法等。


从Python到机器学习

本系列课程将全面介绍深度学习入门的应用知识。包括从Python基础开始，到深度学习框架Tensorflow的使用方法。是一套简练风趣，易懂易学的入门课程。


机器学习思维

本课程围绕机器学习思维，讨论机器学习的运用方法，能力范围，技术种类，以及机器学习与人类学习的不同点。


当复杂网络遇上机器学习

本课程将对网络嵌入、图卷积神经网络、图注意力模型等新技术的算法及应用，以及多种网络嵌入算法，链路预测问题等问题座以一个简单介绍，并对网络上的深度学习问题进行了展望。


本中文词条由Xugami审校，薄荷编辑，如有问题，欢迎在讨论页面留言。


本词条内容源自wikipedia及公开资料，遵守 CC3.0协议。

取自“https://wiki.swarma.org/index.php?title=EM算法&oldid=28024”