初识EM算法

1. 当我们面对“看不见”的数据

在统计建模和机器学习中，我们通常的目标是根据观测到的数据来估计模型的参数。例如，在回归问题中，我们根据样本点 $(x_i,y_i)$ 来估计回归系数 $\beta$；在分类问题中，我们根据带有标签的样本来估计分类边界或概率分布参数。

然而，在许多实际场景中，数据并非总是完整的，或者模型中存在我们无法直接观测到的“内部状态”或“隐藏变量”。

• 数据缺失：比如，问卷调查中有些受访者没有回答所有问题。
• 隐变量：在聚类问题中，每个数据点属于哪个簇，在观测之前是未知的“隐藏”信息。在混合模型（Mixture Model）中，数据是由多个不同的概率分布生成的，但我们并不知道每个观测值具体是由哪个分布生成的。

当数据不完整或存在隐变量时，传统的参数估计方法（如直接应用极大似然估计）会变得非常困难，甚至无法进行。这时，EM（Expectation-Maximization，期望最大化）算法就应运而生，它提供了一种优雅而通用的迭代方法来解决这类问题。

2. 回顾极大似然估计 (Maximum Likelihood Estimation, MLE)

在深入EM算法之前，我们先来复习一下参数估计中的核心方法——极大似然估计（MLE）。

2.1 MLE的直观思想

想象一下，我们有一枚硬币，它可能是不均匀的，正面朝上的概率是 $p$。我们抛了10次硬币，结果是8次正面，2次反面。现在，我们想估计这个 $p$ 值是多少。

• 如果 $p=0.5$，出现8正2反的概率是 $\left(\genfrac{}{}{0ex}{}{10}{8}\right)(0.5{)}^8(0.5{)}^2\approx 0.044$。
• 如果 $p=0.8$，出现8正2反的概率是 $\left(\genfrac{}{}{0ex}{}{10}{8}\right)(0.8{)}^8(0.2{)}^2\approx 0.302$。

显然，在 $p=0.8$ 的情况下，我们观测到8正2反的这个结果的可能性更大。 极大似然估计的核心思想就是：给定了一组已经观测到的数据，我们去寻找这样一组模型参数，使得这组参数下，我们观测到的这些数据发生的概率最大。换句话说，让已经发生的事件在我们的模型参数下显得最合理、最可能发生。

2.2 MLE的数学定义

假设我们有一个数据集 $D=\{x_1,x_2,\dots ,x_N\}$，这些数据是从某个概率分布 $P(x;\theta )$ 中独立同分布（i.i.d.）抽样得到的，其中 $\theta$ 是我们想要估计的未知参数（例如，高斯分布的均值 $\mu$ 和方差 ${\sigma }^2$）。

似然函数（Likelihood Function） $L(\theta ;D)$ 定义为在给定参数 $\theta$ 的情况下，观测到数据集 $D$ 的概率（或概率密度积）：

$L(\theta ;D)=P(x_1,x_2,\dots ,x_N;\theta )$

由于数据是独立同分布的，似然函数可以写成各个数据点概率的乘积：

$L(\theta ;D)={\prod }_{i=1}^{N}P(x_i;\theta )$

我们的目标是找到使似然函数最大的参数 ${\hat{\theta }}_{MLE}$：

${\hat{\theta }}_{MLE}=\arg {\max }_{\theta }L(\theta ;D)=\arg {\max }_{\theta }{\prod }_{i=1}^{N}P(x_i;\theta )$

为了计算方便（将乘积转化为求和，且不改变最大化位置），我们通常最大化对数似然函数（Log-Likelihood Function）：

${\hat{\theta }}_{MLE}=\arg {\max }_{\theta }\log L(\theta ;D)=\arg {\max }_{\theta }{\sum }_{i=1}^N\log P(x_i;\theta )$

求解步骤：通常通过对对数似然函数求导，并令导数等于零来找到参数的估计值。

2.3 完整数据下的MLE示例

我们尝试分析这样一个模型：$X_{ij}=\mu +{\alpha }_i+{\beta }_j+{ϵ}_{ij}$，其中 ${ϵ}_{ij}\sim N(0,{\sigma }^2)$。这是一个常见的两因素方差分析（Two-way ANOVA） 模型或称为 线性模型 。

• $X_{ij}$: 第 $i$ 行第 $j$ 列的观测值。
• $\mu$: 总平均值。
• ${\alpha }_i$: 第 $i$ 行（或因子A的第 $i$ 个水平）的效应。
• ${\beta }_j$: 第 $j$ 列（或因子B的第 $j$ 个水平）的效应。
• ${ϵ}_{ij}$: 误差项，服从均值为0，方差为 ${\sigma }^2$ 的正态分布。

假设所有数据 $X_{ij}$ 都被完整观测到了。我们的目标是估计参数 $\theta =\{\mu ,{\alpha }_1,\dots ,{\alpha }_I,{\beta }_1,\dots ,{\beta }_J,{\sigma }^2\}$。

由于 ${ϵ}_{ij}=X_{ij}-(\mu +{\alpha }_i+{\beta }_j)$，并且 ${ϵ}_{ij}$ 服从正态分布，那么 $X_{ij}$ 也服从正态分布：

$X_{ij}\sim N(\mu +{\alpha }_i+{\beta }_j,{\sigma }^2)$

单个观测值 $X_{ij}$ 的概率密度函数是：

$P(X_{ij};\theta )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(X_{ij}-(\mu +{\alpha }_i+{\beta }_j){)}^2}{2{\sigma }^2}\right)$

对数似然函数 是所有观测值的对数概率密度之和：

$\log L(\theta ;D)={\sum }_{i,j}\log \left(\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(X_{ij}-(\mu +{\alpha }_i+{\beta }_j){)}^2}{2{\sigma }^2}\right)\right)$

$={\sum }_{i,j}\left(-\frac{1}{2}\log (2\pi {\sigma }^2)-\frac{(X_{ij}-(\mu +{\alpha }_i+{\beta }_j){)}^2}{2{\sigma }^2}\right)$

为了最大化这个对数似然函数，我们通常需要最小化误差平方和 ${\sum }_{i,j}(X_{ij}-(\mu +{\alpha }_i+{\beta }_j){)}^2$。而对于这个完整的模型，参数的MLE估计可以通过简单的均值和残差计算得到（例如，$\hat{\mu }={\bar{X}}_{\cdot \cdot }$, ${\hat{\alpha }}_i={\bar{X}}_{i\cdot }-{\bar{X}}_{\cdot \cdot }$, ${\hat{\beta }}_j={\bar{X}}_{\cdot j}-{\bar{X}}_{\cdot \cdot }$）。这些都是解析解，可以直接计算出来。

3. 如果存在不完整数据与隐变量呢？

现在，我们来看当数据不完整或存在隐变量时，极大似然估计会遇到什么困难。

3.1 不完整数据的问题

回到上面的 $X_{ij}=\mu +{\alpha }_i+{\beta }_j+{ϵ}_{ij}$ 模型。 假设 $X_{23}$ =（第二行第三列的观测值）=缺失了，即我们无法观测到它。

• 直接MLE的困难： 由于 $X_{23}$ 未知，我们无法直接计算包含 $X_{23}$ 的对数似然项 $\log P(X_{23};\theta )$。 此时，我们只能对观测到的数据 $D_{obs}$ 进行似然估计。这意味着我们必须把缺失数据 $X_{mis}$ 积分掉：

  $L(\theta ;D_{obs})=P(D_{obs};\theta )=\int P(D_{obs},X_{mis};\theta )d{X}_{mis}$

  这个积分通常非常复杂，甚至没有解析解，导致直接最大化 $L(\theta ;D_{obs})$ 变得异常困难。

3.2 隐变量问题

除了简单的缺失值，还有一类更普遍的问题是模型中包含隐变量（Latent Variables）。这些变量是模型内部的状态，它们影响着观测数据，但我们本身无法直接观测到它们。

经典案例：硬币投掷问题 假设有两枚不均匀的硬币A和硬币B。

• 硬币A：正面朝上的概率是 ${\theta }_A$。
• 硬币B：正面朝上的概率是 ${\theta }_B$。 我们随机选择一枚硬币（但我们不知道选了哪一枚），然后抛掷10次，记录结果。重复这个过程5次。 观测数据 $Y$：5次实验的抛掷结果序列，例如：
• 第一次：HTHTHTHTHH (8次正面，2次反面)
• 第二次：HHHHHHHHHH (10次正面，0次反面)
• … 隐变量 $Z$：每一次实验，我们选择了哪一枚硬币（是A还是B）。这是我们不知道的，但它决定了每次实验结果的生成分布。

我们的目标是根据观测到的抛掷结果，估计两枚硬币的参数 ${\theta }_A$ 和 ${\theta }_B$。

• 直接MLE的困难： 如果每次实验我们都知道选择了哪枚硬币，那估计 ${\theta }_A$ 和 ${\theta }_B$ 就很简单：直接统计硬币A的所有抛掷结果来估计 ${\theta }_A$，统计硬币B的所有抛掷结果来估计 ${\theta }_B$。 但现在，我们不知道。似然函数是关于观测数据的，它需要对隐变量进行求和（或积分）：

  $L(\theta ;Y)=P(Y;\theta )={\sum }_{Z}P(Y,Z;\theta )$

  这个求和（或积分）同样会使优化变得复杂。

3.3 为什么直接MLE失效或困难？

无论是缺失数据还是隐变量，它们都导致我们无法直接得到“完整数据”的似然函数。我们能够写的似然函数是关于观测数据的似然 $P(Y_{obs};\theta )$。

$\log P(Y_{obs};\theta )=\log {\sum }_{Z}P(Y_{obs},Z;\theta )$

注意 $\log {\sum }_Z\dots$ 中的求和是在对数内部。由于对数的非线性，这个求和使得对数似然函数通常是非凸的，并且没有解析解。对其求导并置零不再简单，甚至无法直接应用。 EM算法正是为了解决这类问题而设计的。

对于这两个问题，在 【【合集】十分钟 机器学习 系列视频 《统计学习方法》——监督学习篇】 中详细的讲解了两个例子，有兴趣可以去看看。

4. 概率论基础

为了理解EM算法，我们需要回顾一些核心的概率论概念。

4.1 基本概念：联合概率与条件概率

• 联合概率（Joint Probability）：$P(A,B)$。表示事件A和事件B同时发生的概率。 例如：$P(\text{硬币A},\text{正面})$ 表示选择硬币A并且抛出正面的概率。
• 条件概率（Conditional Probability）：$P(A|B)$。表示在事件B已经发生的条件下，事件A发生的概率。 定义式：$P(A|B)=\frac{P(A,B)}{P(B)}$，前提是 $P(B)>0$。 由定义可得乘法公式：$P(A,B)=P(A|B)P(B)=P(B|A)P(A)$。 例如：$P(\text{正面}|\text{硬币A})$ 表示在选择了硬币A的条件下，抛出正面的概率，这个就是硬币A的偏置 ${\theta }_A$。

4.2 贝叶斯公式 (Bayes’ Formula)

贝叶斯公式是概率论中一个非常重要的公式，它描述了在已知某些先验信息和观测数据的情况下，如何更新我们对事件发生的信念（后验概率）。

公式形式：

$P(A|B)=\frac{P(B|A)P(A)}{P(B)}$

其中：

• $P(A|B)$ 是后验概率（Posterior Probability）：在观察到事件B后，事件A发生的概率。这是我们最关心的，因为我们想根据观测数据来更新对参数或隐变量的认识。
• $P(B|A)$ 是似然（Likelihood）：在事件A发生的条件下，观察到事件B的概率。
• $P(A)$ 是先验概率（Prior Probability）：在观察到事件B之前，事件A发生的概率。
• $P(B)$ 是边际似然（Marginal Likelihood） 或证据（Evidence）：事件B发生的总概率。通常通过对所有可能的A进行求和得到：$P(B)={\sum }_{A}P(B|A)P(A)$。

直观理解贝叶斯公式： 贝叶斯公式提供了一个框架，让我们能够根据新的证据（观测数据B） 来更新我们对某个假设（事件A） 的信念 。

• 我们有一个关于A的初始信念（先验概率 $P(A)$）。
• 我们观察到了一些数据B。
• 贝叶斯公式告诉我们，如何将 $P(A)$ 结合数据B的证据 $P(B|A)$，来得到一个更精确的信念（后验概率 $P(A|B)$）。

在EM算法的E步中，我们将大量使用贝叶斯公式来计算隐变量的后验概率。

5. EM算法：期望最大化 (Expectation-Maximization)

5.1 EM算法的迭代“猜测”与“优化”

EM算法是为了解决包含隐变量或不完整数据的MLE问题而提出的。它是一个迭代过程，在每一步中，它都包含两个子步骤：期望步（E-step） 和最大化步（M-step）。

想象我们正在玩一个拼图游戏，但有一些拼图块被遮住了（隐变量）。

1. E步（Expectation）—— 猜测阶段： 我们根据目前对整幅图的猜测（模型参数），去猜测那些被遮住的拼图块可能是什么样子，或者它们有多大的概率是某种样子。我们不是确定地填上它们，而是给出它们所有可能的“期望”形式。 （例如：根据当前估计的硬币偏置，计算每次实验是硬币A还是硬币B抛出的概率。）
2. M步（Maximization）—— 优化阶段： 现在，我们拥有了对所有拼图块的“完整”信息（包括我们猜测的那些）。基于这些“完整”信息，我们去优化我们的整幅图的猜测，使得它最符合我们现在所看到的所有信息（包括被遮住的那些的“期望”）。 （例如：根据E步计算出的每次实验是A还是B的“概率”，重新估计硬币A和硬币B的最佳偏置。）

这个“猜测”和“优化”的过程会交替进行，每次迭代都会使模型的似然函数单调增加，直到收敛。

5.2 EM算法的正式步骤

EM算法的输入通常是观察变量数据 $Y$、隐变量数据 $Z$ 的结构，以及联合概率分布 $P(Y,Z|\theta )$（也称为完整数据概率），输出是模型参数 $\theta$ 的估计值。

算法 9.1 EM算法

输入： 观测变量数据 $Y$，隐变量数据 $Z$，联合概率分布 $P(Y,Z|\theta )$，条件分布 $P(Z|Y,\theta )$。 输出： 模型参数 $\theta$ 的估计值。

(1) 选择参数的初始值 ${\theta }^{(0)}$，开始迭代。

(2) E步（Expectation Step）：记 ${\theta }^{(t)}$ 为第 $t$ 次迭代的参数估计值。在第 $t+1$ 次迭代的E步中，计算 Q 函数：

$Q(\theta ,{\theta }^{(t)})=E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(t)}]$ $={\sum }_Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(t)})(9.9)$

• 解释：

  • $Q(\theta ,{\theta }^{(t)})$ 是一个关于新参数 $\theta$ 的函数，而 ${\theta }^{(t)}$ 是固定的旧参数。

  • $P(Y,Z|\theta )$ 是完整数据 的对数似然函数。如果 $Z$ 是已知的，我们就可以直接最大化这个函数。

  • $E_Z[\dots |Y,{\theta }^{(t)}]$ 表示对隐变量 $Z$ 的期望。这个期望是在给定观测数据 $Y$ 和当前参数估计 ${\theta }^{(t)}$ 的条件下进行的。

  • 具体地，$P(Z|Y,{\theta }^{(t)})$ 是隐变量 $Z$ 的后验概率，它通过贝叶斯公式计算得到： $P(Z|Y,{\theta }^{(t)})=\frac{P(Y,Z|{\theta }^{(t)})}{{\sum }_{Z}P(Y,Z|{\theta }^{(t)})}=\frac{P(Y|Z,{\theta }^{(t)})P(Z|{\theta }^{(t)})}{P(Y|{\theta }^{(t)})}$

    这个后验概率 $P(Z|Y,{\theta }^{(t)})$ 就是E步中对隐变量 $Z$ 的“猜测”或“概率分布”。

(3) M步（Maximization Step）：求使 $Q(\theta ,{\theta }^{(t)})$ 极大化的 $\theta$，确定第 $t+1$ 次迭代的参数估计值 ${\theta }^{(t+1)}$。

${\theta }^{(t+1)}=\arg {\max }_{\theta }Q(\theta ,{\theta }^{(t)})(9.10)$

• 解释：
  • 在E步中，我们已经计算出了Q函数，它实际上是一个关于 $\theta$ 的函数，并且其中的隐变量 $Z$ 已经被其期望代替了。
  • M步就是简单地对这个Q函数进行最大化，就像进行一次标准的极大似然估计一样。由于Q函数通常比原始的观测数据对数似然函数更容易最大化，所以这一步是可行的。这通常涉及到求导并令导数等于零。

(4) 重复步骤 (2) 和步骤 (3)，直到收敛。

• 收敛准则：当参数估计值 ${\theta }^{(t+1)}$ 和 ${\theta }^{(t)}$ 之间的差异非常小，或者对数似然函数 $\log P(Y|{\theta }^{(t)})$ 的增长非常小，就可以认为算法收敛。

5.3 E步 (Expectation Step)：计算期望的完整数据对数似然（Q函数）

E步的核心是计算Q函数。Q函数是完整数据对数似然函数 $\log P(Y,Z|\theta )$ 关于隐变量 $Z$ 的期望。这个期望是在给定观测数据 $Y$ 和当前参数 ${\theta }^{(t)}$ 的条件下计算的。

• 为什么需要期望？ 因为隐变量 $Z$ 是未知的，我们无法直接计算 $P(Y,Z|\theta )$。因此，我们用它的期望值来代替。这个期望值是基于我们对隐变量最合理的“猜测”——即它们的后验概率分布 $P(Z|Y,{\theta }^{(t)})$。
• Q函数的意义：Q函数可以看作是“填补了缺失信息”的对数似然函数。我们用 $P(Z|Y,{\theta }^{(t)})$ 作为权重，对所有可能的 $Z$ 值对应的 $\log P(Y,Z|\theta )$ 进行加权求和（对于离散隐变量）或积分（对于连续隐变量）。

5.4 M步 (Maximization Step)：最大化Q函数

M步的核心是最大化E步得到的Q函数。这一步的目标是找到新的参数 ${\theta }^{(t+1)}$，使得Q函数取最大值。

• 为什么可以最大化？ 在E步中，Q函数已经将隐变量“消除”了（通过求期望），使得它成为一个关于 $\theta$ 的纯函数。通常，这个函数比原始的观测数据对数似然函数更容易最大化，甚至常常有解析解。
• 与MLE的关系：M步本质上就是在进行一次标准的极大似然估计，但它是在一个“填补了缺失信息”（或期望了隐变量）的对数似然函数上进行的。

5.5 EM算法的收敛性：为什么它能工作？

EM算法的一个重要性质是，它能保证每一步迭代都使观测数据的对数似然函数 $P(Y|\theta )$ 单调不减。这意味着，算法最终会收敛到局部最优解。

• Jensen’s Inequality（琴生不等式）的作用： 它是EM算法收敛性证明的核心工具。 Jensen不等式：对于一个凹函数 $f$（例如 $\log (x)$），有 $E[f(X)]\le f(E[X])$。对于凸函数 $f$，有 $E[f(X)]\ge f(E[X])$。 在EM算法的证明中，我们会利用 $\log (x)$ 是凹函数这个性质，推导出 $L(\theta )$ 会单调递增。因此EM算法的每次迭代都在不断提升观测数据的似然值。它不会让模型的性能变差，只会让它变得更好或保持不变，直到达到一个稳定点（局部最优解）。

6. EM算法的经典案例——硬币投掷问题

问题设定：

• 有两枚硬币A和B，它们的正面朝上概率分别为 ${\theta }_A$ 和 ${\theta }_B$。
• 我们进行5次独立实验。每次实验先随机选择一枚硬币（假设选择A或B的概率都是0.5，这是我们模型的一部分），然后抛掷10次。
• 观测数据 $Y$：5次实验的抛掷结果，用 $y_j$ 表示第 $j$ 次实验的正面次数和反面次数。 例如，第1次实验：$(H=8,T=2)$；第2次实验：$(H=10,T=0)$；…；第5次实验：$(H=5,T=5)$。
• 隐变量 $Z$：第 $j$ 次实验选择了哪枚硬币（$Z_j\in \{\text{A, B}\}$）。我们不知道 $Z_j$。
• 目标：估计 ${\theta }_A$ 和 ${\theta }_B$。

完整数据 log-likelihood： 如果我们知道每次选择了哪枚硬币 $Z_j$，那么完整数据的对数似然函数为：

$\log P(Y,Z|{\theta }_A,{\theta }_B)={\sum }_{j=1}^5\left(\log P(Z_j|\theta )+\log P(y_j|Z_j,\theta )\right)$

其中 $P(Z_j|\theta )$ 是选择硬币 $Z_j$ 的概率 (假设为0.5)。 $\log P(y_j|Z_j,\theta )$ 是在给定硬币 $Z_j$ 和参数 $\theta$ 下观测到 $y_j$ 的概率。 例如，如果第 $j$ 次实验选择了硬币A，$y_j$ 有 $k_j$ 个正面和 $10-k_j$ 个反面，那么 $\log P(y_j|\text{A},{\theta }_A)=\log (\genfrac{}{}{0ex}{}{10}{k_j})+k_j\log {\theta }_A+(10-k_j)\log (1-{\theta }_A)$。

6.1 EM算法的迭代计算演示 (第一轮)

步骤 0：初始化参数 ${\theta }^{(0)}$ 我们随机猜测两枚硬币的初始偏置： ${\theta }_A^{(0)}=0.60$ ${\theta }_B^{(0)}=0.50$

步骤 1：E步 (Expectation Step)

核心：计算在当前参数 ${\theta }_A^{(0)},{\theta }_B^{(0)}$ 下，每次实验是由硬币A还是硬币B抛出的后验概率 $P(Z_j|y_j,{\theta }^{(0)})$。

假设第 $j$ 次实验观测到 $k_j$ 次正面（H）和 $10-k_j$ 次反面（T）。 我们使用贝叶斯公式计算 $P(\text{A}|y_j,{\theta }^{(0)})$ 和 $P(\text{B}|y_j,{\theta }^{(0)})$。

$P(\text{A}|y_j,{\theta }^{(0)})=\frac{P(y_j|\text{A},{\theta }_A^{(0)})P(\text{A})}{P(y_j|\text{A},{\theta }_A^{(0)})P(\text{A})+P(y_j|\text{B},{\theta }_B^{(0)})P(\text{B})}$

由于 $P(\text{A})=P(\text{B})=0.5$（选择硬币A或B的先验概率），公式简化为：

$P(\text{A}|y_j,{\theta }^{(0)})=\frac{P(y_j|\text{A},{\theta }_A^{(0)})}{P(y_j|\text{A},{\theta }_A^{(0)})+P(y_j|\text{B},{\theta }_B^{(0)})}$

其中 $P(y_j|\text{Coin},{\theta }_{\text{Coin}})=(\genfrac{}{}{0ex}{}{10}{k_j}){\theta }_{\text{Coin}}^{k_j}(1-{\theta }_{\text{Coin}}{)}^{10-k_j}$。

现在 让我们计算的第一个例子：第1次实验，8H, 2T。

• 在 ${\theta }_A^{(0)}=0.60$ 下，得到8H2T的似然：$P(y_1|\text{A},{\theta }_A^{(0)})=(\genfrac{}{}{0ex}{}{10}{8})(0.60{)}^8(0.40{)}^2\approx 0.0450$
• 在 ${\theta }_B^{(0)}=0.50$ 下，得到8H2T的似然：$P(y_1|\text{B},{\theta }_B^{(0)})=(\genfrac{}{}{0ex}{}{10}{8})(0.50{)}^8(0.50{)}^2\approx 0.0439$

计算第1次实验是硬币A抛出的后验概率（称为“责任”或“responsibility”）： $P(\text{A}|y_1,{\theta }^{(0)})=\frac{0.0450}{0.0450+0.0439}\approx \frac{0.0450}{0.0889}\approx 0.506$ 同理，第1次实验是硬币B抛出的后验概率： $P(\text{B}|y_1,{\theta }^{(0)})=1-0.506=0.494$ 这表示第一次实验有约50.6%的可能是由硬币A抛出的。

对所有5次实验重复E步： 例如，表格中“Coin A”列的 $2.2H,2.8T$ 并不是实际的观测值，而是对硬币A的总贡献，即该实验被认为是硬币A抛出的期望正面/反面次数。

• 第1次实验 (8H, 2T): 认为是A的责任是0.506
  • A的期望贡献：$0.506\times 8\text{H}=4.048\text{H}$，$0.506\times 2\text{T}=1.012\text{T}$
• 第2次实验 (10H, 0T): 责任是0.80
  • A的期望贡献：$0.80\times 10\text{H}=8.0\text{H}$，$0.80\times 0\text{T}=0\text{T}$
• 第3次实验 (7H, 3T): 责任是0.73
  • A的期望贡献：$0.73\times 7\text{H}=5.11\text{H}$，$0.73\times 3\text{T}=2.19\text{T}$
• 第4次实验 (3H, 7T): 责任是0.35
  • A的期望贡献：$0.35\times 3\text{H}=1.05\text{H}$，$0.35\times 7\text{T}=2.45\text{T}$
• 第5次实验 (5H, 5T): 责任是0.65
  • A的期望贡献：$0.65\times 5\text{H}=3.25\text{H}$，$0.65\times 5\text{T}=3.25\text{T}$

将所有实验中硬币A的期望贡献加起来：

• A的总期望正面数：$4.048+8.0+5.11+1.05+3.25=21.458\approx 21.3$
• A的总期望反面数：$1.012+0+2.19+2.45+3.25=8.902\approx 8.7$

同样对硬币B进行计算：

• 第1次实验 (8H, 2T): 认为是B的责任是0.494
  • B的期望贡献：$0.494\times 8\text{H}=3.952\text{H}$，$0.494\times 2\text{T}=0.988\text{T}$ … (计算所有5次实验对B的期望贡献)
• B的总期望正面数：$3.952+\cdots \approx 11.7$
• B的总期望反面数：$0.988+\cdots \approx 8.4$

这些期望的正面/反面计数就是E步的输出，它们构成了Q函数的一部分。

步骤 2：M步 (Maximization Step)

核心：使用E步中计算出的“期望的完整数据”，重新估计参数 ${\theta }_A$ 和 ${\theta }_B$。 这就像是，我们现在知道了（以概率形式）每一枚硬币贡献了多少次正面和反面，于是我们就可以像在完整数据下一样，简单地用频率来估计概率。

更新 ${\theta }_A^{(1)}$：用所有实验中硬币A的总期望正面数除以硬币A的总期望抛掷次数（正面+反面）。

${\theta }_A^{(1)}=\frac{\text{A的总期望正面数}}{\text{A的总期望正面数}+\text{A的总期望反面数}}=\frac{21.3}{21.3+8.7}=\frac{21.3}{30}\approx 0.71$

更新 ${\theta }_B^{(1)}$：

${\theta }_B^{(1)}=\frac{\text{B的总期望正面数}}{\text{B的总期望正面数}+\text{B的总期望反面数}}=\frac{11.7}{11.7+8.4}=\frac{11.7}{20.1}\approx 0.58$

6.2 迭代过程与收敛

将新的参数估计 ${\theta }_A^{(1)}=0.71,{\theta }_B^{(1)}=0.58$ 作为下一次迭代的初始值 ${\theta }^{(1)}$，然后重复E步和M步。

第二轮E步：使用 ${\theta }_A^{(1)}=0.71,{\theta }_B^{(1)}=0.58$ 重新计算每次实验是硬币A还是硬币B抛出的后验概率。这些概率会发生变化，因为我们对硬币偏置的猜测更接近真实了。 第二轮M步：使用这些新的后验概率，再次计算A和B的总期望正面/反面数，并更新 ${\theta }_A^{(2)}$ 和 ${\theta }_B^{(2)}$。

这个过程会一直重复，直到参数 ${\theta }_A$ 和 ${\theta }_B$ 的估计值趋于稳定，不再发生显著变化。算法最终会收敛到一个局部最优解。

好的，我们继续深入学习EM算法，并以经典的三硬币模型为例，详细解析其在实际问题中的应用。这个例子非常典型，能够帮助您巩固对EM算法E步和M步的理解。

---

7. 三硬币模型实例分析

7.1 模型设定与问题描述

问题描述： 假设有三枚硬币，分别记作A、B、C。

• 硬币A：正面朝上的概率为 $\pi$。
• 硬币B：正面朝上的概率为 $p$。
• 硬币C：正面朝上的概率为 $q$。

实验过程： 我们进行 $N$ 次（假设 $N=10$）独立的重复试验。每次试验的步骤如下：

1. 先掷硬币A：根据硬币A的结果决定接下来掷哪枚硬币。
   • 如果硬币A为正面（概率为 $\pi$），则选择硬币B。
   • 如果硬币A为反面（概率为 $1-\pi$），则选择硬币C。
2. 掷选定的硬币：掷选定的硬币（B或C）1次，并记录其结果（正面为1，反面为0）。

观测数据 $Y$：我们只能观测到最终掷出硬币（B或C）的结果。 例如， $Y=(1,1,0,1,0,0,1,0,1,1)$，表示10次试验中最终观测到的结果。 隐变量 $Z$：每次试验中，我们掷的是硬币B还是硬币C，这是我们无法观测到的。

目标：根据观测数据 $Y$，估计三硬币模型的参数 $\theta =(\pi ,p,q)$。

7.2 模型结构与参数定义

• 观测变量 $Y$：每次试验的最终结果 $y_j\in \{0,1\}$，$j=1,\dots ,N$。
• 隐变量 $Z$：每次试验中选择的硬币 $z_j\in \{B,C\}$，$j=1,\dots ,N$。
• 模型参数 $\theta =(\pi ,p,q)$。

7.3 完整数据下的概率分布

如果我们能够观测到隐变量 $Z$（即知道每次试验掷的是B还是C），那么完整数据为 $(Y,Z)=\{(y_j,z_j){\}}_{j=1}^N$。

对于单次试验 $j$，完整数据的联合概率 $P(y_j,z_j|\theta )$ 可以通过乘法公式得到：

$P(y_j,z_j|\theta )=P(y_j|z_j,\theta )P(z_j|\theta )$

• $P(z_j|\theta )$：这是根据硬币A的结果决定选择硬币B或C的概率。
  • 如果 $z_j=B$，则 $P(z_j=B|\theta )=\pi$。
  • 如果 $z_j=C$，则 $P(z_j=C|\theta )=1-\pi$。
• $P(y_j|z_j,\theta )$：在选择了硬币 $z_j$ 的条件下，掷出 $y_j$ 的概率。
  • 如果 $z_j=B$：$P(y_j|z_j=B,\theta )=p^{y_j}(1-p{)}^{1-y_j}$。（伯努利分布的PMF：若 $y_j=1$ 则为 $p$，若 $y_j=0$ 则为 $1-p$）
  • 如果 $z_j=C$：$P(y_j|z_j=C,\theta )=q^{y_j}(1-q{)}^{1-y_j}$。

所以，完整数据的联合概率可以写为：

• 当 $z_j=B$ 时：$P(y_j,z_j=B|\theta )=p^{y_j}(1-p{)}^{1-y_j}\cdot \pi$
• 当 $z_j=C$ 时：$P(y_j,z_j=C|\theta )=q^{y_j}(1-q{)}^{1-y_j}\cdot (1-\pi )$

完整数据的对数似然函数是所有试验的对数联合概率之和：

$\log P(Y,Z|\theta )={\sum }_{j=1}^N\log P(y_j,z_j|\theta )$ $={\sum }_{j=1}^N[z_j\log (\pi p^{y_j}(1-p{)}^{1-y_j})+(1-z_j)\log ((1-\pi )q^{y_j}(1-q{)}^{1-y_j})]$

这里我们用指示变量 $z_j$ 来表示隐变量，如果 $z_j=B$ 设为1，如果 $z_j=C$ 设为0。

7.4 EM算法步骤详解

我们将按照EM算法的E步和M步进行迭代。

步骤 0：初始化参数 随机选择参数的初始值 ${\theta }^{(0)}=({\pi }^{(0)},p^{(0)},q^{(0)})$。

步骤 1：E步

核心任务：计算在给定观测数据 $y_j$ 和当前参数 ${\theta }^{(t)}$ 的条件下，隐变量 $z_j$ 的后验概率 $P(z_j|y_j,{\theta }^{(t)})$。 对于每一次试验 $j$，我们计算它是由硬币B抛出的概率，记为 ${\mu }_j^{(t+1)}$：

${\mu }_j^{(t+1)}=P(z_j=B|y_j,{\theta }^{(t)})$

根据贝叶斯公式，我们有：

$P(z_j=B|y_j,{\theta }^{(t)})=\frac{P(y_j|z_j=B,{\theta }^{(t)})P(z_j=B|{\theta }^{(t)})}{P(y_j|z_j=B,{\theta }^{(t)})P(z_j=B|{\theta }^{(t)})+P(y_j|z_j=C,{\theta }^{(t)})P(z_j=C|{\theta }^{(t)})}$

展开各项：

• $P(y_j|z_j=B,{\theta }^{(t)})=(p^{(t)}{)}^{y_j}(1-p^{(t)}{)}^{1-y_j}$
• $P(z_j=B|{\theta }^{(t)})={\pi }^{(t)}$
• $P(y_j|z_j=C,{\theta }^{(t)})=(q^{(t)}{)}^{y_j}(1-q^{(t)}{)}^{1-y_j}$
• $P(z_j=C|{\theta }^{(t)})=1-{\pi }^{(t)}$

代入得到：

${\mu }_j^{(t+1)}=\frac{(p^{(t)}{)}^{y_j}(1-p^{(t)}{)}^{1-y_j}\cdot {\pi }^{(t)}}{(p^{(t)}{)}^{y_j}(1-p^{(t)}{)}^{1-y_j}\cdot {\pi }^{(t)}+(q^{(t)}{)}^{y_j}(1-q^{(t)}{)}^{1-y_j}\cdot (1-{\pi }^{(t)})}$

• 解释：这个 ${\mu }_j^{(t+1)}$ 代表了在第 $t$ 轮参数估计下，第 $j$ 次观测结果 $y_j$ 是由硬币B（而不是C）产生的“责任”或“概率”。这个值就是我们对隐变量 $z_j$ 的“猜测”。

同时，我们可以得到由硬币C抛出的概率 $P(z_j=C|y_j,{\theta }^{(t)})=1-{\mu }_j^{(t+1)}$。

Q函数的构建： Q函数是完整数据对数似然函数 $\log P(Y,Z|\theta )$ 关于隐变量 $Z$ 的期望。 将 $Z$ 的值替换为其后验概率 ${\mu }_j^{(t+1)}$：

$$\begin{array}{rl}Q(\theta ,{\theta }^{(t)})& =\sum _{j=1}^N{E}_{z_j}[\log P(y_j,z_j|\theta )|y_j,{\theta }^{(t)}]\\ & =\sum _{j=1}^N\left[{\mu }_j^{(t+1)}\log P(y_j,z_j=B|\theta )+(1-{\mu }_j^{(t+1)})\log P(y_j,z_j=C|\theta )\right]\end{array}$$

将具体的概率分布代入：

$$\begin{array}{rl}Q(\theta ,{\theta }^{(t)})=& \sum _{j=1}^N[{\mu }_j^{(t+1)}\left(\log \pi +y_j\log p+(1-y_j)\log (1-p)\right)\\ & +(1-{\mu }_j^{(t+1)})\left(\log (1-\pi )+y_j\log q+(1-y_j)\log (1-q)\right)]\end{array}$$

• 解释：Q函数将我们无法观测的隐变量 $z_j$ 用其期望的后验概率 ${\mu }_j^{(t+1)}$ 进行了加权。现在Q函数完全是关于参数 $\theta =(\pi ,p,q)$ 的显式函数。

步骤 2：M步（Maximization Step）

核心任务：最大化Q函数 $Q(\theta ,{\theta }^{(t)})$ 关于新参数 $\theta =(\pi ,p,q)$。 这一步通过对Q函数求偏导并令其为0来完成。我们将Q函数分别对 $\pi ,p,q$ 求偏导。

1. 估计 ${\pi }^{(t+1)}$： 只看Q函数中包含 $\pi$ 的项：

${\sum }_{j=1}^N\left[{\mu }_j^{(t+1)}\log \pi +(1-{\mu }_j^{(t+1)})\log (1-\pi )\right]$

对 $\pi$ 求偏导并令其为0：

$$\begin{array}{rl}\frac{\partial Q}{\partial \pi }& =\sum _{j=1}^N\left[\frac{{\mu }_j^{(t+1)}}{\pi }-\frac{1-{\mu }_j^{(t+1)}}{1-\pi }\right]=0\\ \frac{1}{\pi }\sum _{j=1}^N{\mu }_j^{(t+1)}& =\frac{1}{1-\pi }\sum _{j=1}^N(1-{\mu }_j^{(t+1)})\\ (1-\pi )\sum _{j=1}^N{\mu }_j^{(t+1)}& =\pi \left(N-\sum _{j=1}^N{\mu }_j^{(t+1)}\right)\\ \sum _{j=1}^N{\mu }_j^{(t+1)}-\pi \sum _{j=1}^N{\mu }_j^{(t+1)}& =N\pi -\pi \sum _{j=1}^N{\mu }_j^{(t+1)}\\ \sum _{j=1}^N{\mu }_j^{(t+1)}& =N\pi \end{array}$$

所以，$\pi$ 的更新公式为：

${\pi }^{(t+1)}=\frac{{\sum }_{j=1}^N{\mu }_j^{(t+1)}}{N}$

• 新估计的 $\pi$ 等于所有试验中由硬币B抛出的期望次数除以总试验次数。这符合MLE的直觉：用频率近似概率。

2. 估计 $p^{(t+1)}$： 只看Q函数中包含 $p$ 的项： ${\sum }_{j=1}^N{\mu }_j^{(t+1)}\left(y_j\log p+(1-y_j)\log (1-p)\right)$ 对 $p$ 求偏导并令其为0：

$$\begin{array}{rl}\frac{\partial Q}{\partial p}& =\sum _{j=1}^N{\mu }_j^{(t+1)}\left(\frac{y_j}{p}-\frac{1-y_j}{1-p}\right)=0\\ \frac{1}{p}\sum _{j=1}^N{\mu }_j^{(t+1)}{y}_j& =\frac{1}{1-p}\sum _{j=1}^N{\mu }_j^{(t+1)}(1-y_j)\\ (1-p)\sum _{j=1}^N{\mu }_j^{(t+1)}{y}_j& =p\sum _{j=1}^N{\mu }_j^{(t+1)}-p\sum _{j=1}^N{\mu }_j^{(t+1)}{y}_j\\ \sum _{j=1}^N{\mu }_j^{(t+1)}{y}_j-p\sum _{j=1}^N{\mu }_j^{(t+1)}{y}_j& =p\sum _{j=1}^N{\mu }_j^{(t+1)}-p\sum _{j=1}^N{\mu }_j^{(t+1)}{y}_j\\ \sum _{j=1}^N{\mu }_j^{(t+1)}{y}_j& =p\sum _{j=1}^N{\mu }_j^{(t+1)}\end{array}$$

所以，$p$ 的更新公式为：

$p^{(t+1)}=\frac{{\sum }_{j=1}^N{\mu }_j^{(t+1)}{y}_j}{{\sum }_{j=1}^N{\mu }_j^{(t+1)}}$

• 新估计的 $p$ 等于所有由硬币B抛出的期望正面次数除以所有由硬币B抛出的期望总抛掷次数。这也是 MLE 的加权频率形式。

3. 估计 $q^{(t+1)}$： 同理，对 $q$ 求偏导并令其为0，得到：

$q^{(t+1)}=\frac{{\sum }_{j=1}^N(1-{\mu }_j^{(t+1)})y_j}{{\sum }_{j=1}^N(1-{\mu }_j^{(t+1)})}$

• 新估计的 $q$ 等于所有由硬币C抛出的期望正面次数除以所有由硬币C抛出的期望总抛掷次数。

步骤 3：重复迭代 使用新的参数估计 ${\theta }^{(t+1)}=({\pi }^{(t+1)},p^{(t+1)},q^{(t+1)})$ 作为下一轮的 ${\theta }^{(t)}$，然后重复E步和M步，直到参数收敛。

7.5 实例数值计算（第一轮迭代）

我们假设观测数据 $Y=(1,1,0,1,0,0,1,0,1,1)$，$N=10$。 初始参数 ${\theta }^{(0)}$：

• ${\pi }^{(0)}=0.5$
• $p^{(0)}=0.5$
• $q^{(0)}=0.5$

第一轮 E步：计算 ${\mu }_j^{(1)}$

对于每个 $j=1,\dots ,10$，计算 ${\mu }_j^{(1)}$：

${\mu }_j^{(1)}=\frac{(p^{(0)}{)}^{y_j}(1-p^{(0)}{)}^{1-y_j}\cdot {\pi }^{(0)}}{(p^{(0)}{)}^{y_j}(1-p^{(0)}{)}^{1-y_j}\cdot {\pi }^{(0)}+(q^{(0)}{)}^{y_j}(1-q^{(0)}{)}^{1-y_j}\cdot (1-{\pi }^{(0)})}$

由于初始参数都是0.5，所以：

• $p^{(0)}=0.5,1-p^{(0)}=0.5$
• $q^{(0)}=0.5,1-q^{(0)}=0.5$
• ${\pi }^{(0)}=0.5,1-{\pi }^{(0)}=0.5$

因此，分子中的 $(p^{(0)}{)}^{y_j}(1-p^{(0)}{)}^{1-y_j}\cdot {\pi }^{(0)}$ 简化为 $(0.5{)}^{y_j}(0.5{)}^{1-y_j}\cdot 0.5=0.5\cdot 0.5=0.25$。 分母中的 $(q^{(0)}{)}^{y_j}(1-q^{(0)}{)}^{1-y_j}\cdot (1-{\pi }^{(0)})$ 同样简化为 $0.5\cdot 0.5=0.25$。

所以，在第一轮E步，所有 ${\mu }_j^{(1)}$ 都是：

${\mu }_j^{(1)}=\frac{0.25}{0.25+0.25}=\frac{0.25}{0.5}=0.5$

• 由于初始参数都是0.5，模型认为每次试验选择B或C的概率是相等的，且B和C抛出正反面的概率也相等。因此，在观测到结果之前，模型认为每次试验由B或C抛出的概率都是0.5。

第一轮 M步：更新参数 ${\theta }^{(1)}$

• 更新 ${\pi }^{(1)}$： ${\pi }^{(1)}=\frac{{\sum }_{j=1}^{10}{\mu }_j^{(1)}}{10}=\frac{10\times 0.5}{10}=0.5$

• 更新 $p^{(1)}$： 观测数据 $Y=(1,1,0,1,0,0,1,0,1,1)$。 其中 $y_j=1$ 的有6次， $y_j=0$ 的有4次。

  ${\sum }_{j=1}^{10}{\mu }_j^{(1)}{y}_j=(0.5\times 1)+(0.5\times 1)+(0.5\times 0)+\cdots +(0.5\times 1)=0.5\times (\text{1的个数})=0.5\times 6=3$

  ${\sum }_{j=1}^{10}{\mu }_j^{(1)}=0.5\times 10=5$

  $p^{(1)}=\frac{3}{5}=0.6$

• 更新 $q^{(1)}$：

  ${\sum }_{j=1}^{10}(1-{\mu }_j^{(1)})y_j=(0.5\times 1)+(0.5\times 1)+\cdots =0.5\times 6=3$

  ${\sum }_{j=1}^{10}(1-{\mu }_j^{(1)})=0.5\times 10=5$

  $q^{(1)}=\frac{3}{5}=0.6$

所以，第一轮迭代后，参数更新为 ${\theta }^{(1)}=(0.5,0.6,0.6)$。

继续迭代：

• 第二轮 E步：现在使用新的参数 ${\theta }^{(1)}=(0.5,0.6,0.6)$ 来计算新的 ${\mu }_j^{(2)}$。
  • 此时，由于 $p^{(1)}\ne q^{(1)}$， ${\mu }_j^{(2)}$ 将不再都是0.5。它们会根据 $y_j$ 是1还是0而有所不同。
    • 如果 $y_j=1$（正面），硬币B和C都可能抛出，但因为 $p^{(1)}=0.6>q^{(1)}=0.6$ (实际上这里相等)，所以 $P(y_j=1|B)$ 和 $P(y_j=1|C)$ 在这一轮相等。这将导致 ${\mu }_j^{(2)}$ 仍然是0.5。
    • 注意：在实际的三硬币模型教学例子中，$p$ 和 $q$ 的初始值通常设置为不同，例如 $p^{(0)}=0.6,q^{(0)}=0.7$，这样在第一轮迭代后，${\mu }_j^{(1)}$ 就会开始分化。

这个迭代过程会继续进行，直到参数收敛到一个稳定值。