期望最大化(EM)算法在伯努利分布混合模型(Mixture of Bernoulli)中的参数估计过程
题目描述:
假设我们有一组高维二值向量(例如,由0和1组成的多维度特征向量,可类比为多张黑白图片的二值化像素),并且我们认为这些数据来自K个不同的伯努利分布混合而成的生成模型。也就是说,每个观测向量是由某个潜在的、不可观测的类别(或组件)生成的,而每个类别都对应一个独立的伯努利分布参数向量。我们的目标是:在给定观测数据但不知道每个数据点属于哪个类别(隐变量)的情况下,估计出这个伯努利混合模型中所有伯努利分布的参数以及每个类别的混合权重。这个问题可以用期望最大化(EM)算法来解决。我们将逐步推导EM算法在此模型中的具体计算公式和迭代过程。
解题过程:
我们将循序渐进地推导,分为五个主要步骤:1) 模型设定与符号定义;2) 完全数据的对数似然;3) E步:计算隐变量的后验期望(责任值);4) M步:最大化期望似然,更新模型参数;5) 算法流程总结。
1. 模型设定与符号定义
- 设有 \(N\) 个观测数据点 \(\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, ..., \mathbf{x}_N\}\),其中每个 \(\mathbf{x}_n\) 是一个 \(D\) 维的二值向量,即 \(\mathbf{x}_n = (x_{n1}, x_{n2}, ..., x_{nD})\),且 \(x_{nd} \in \{0, 1\}\)。
- 假设数据由 \(K\) 个伯努利分布混合生成。我们引入一个隐变量 \(\mathbf{z}_n = (z_{n1}, z_{n2}, ..., z_{nK})\),它是一个 \(K\) 维的one-hot向量,表示数据点 \(n\) 来自哪个混合组件。即,如果 \(\mathbf{x}_n\) 来自第 \(k\) 个组件,则 \(z_{nk} = 1\),其余为0。
- 模型参数包括:
- 混合权重 \(\boldsymbol{\pi} = (\pi_1, \pi_2, ..., \pi_K)\),满足 \(\pi_k \ge 0\) 且 \(\sum_{k=1}^K \pi_k = 1\)。它表示选择第 \(k\) 个组件的先验概率:\(p(z_{nk} = 1) = \pi_k\)。
- 每个组件 \(k\) 对应一个伯努利参数向量 \(\boldsymbol{\mu}_k = (\mu_{k1}, \mu_{k2}, ..., \mu_{kD})\),其中 \(\mu_{kd} \in (0, 1)\) 表示在第 \(k\) 个组件下,第 \(d\) 维特征取值为1的概率。
- 因此,模型的生成过程为:
- 根据 \(\boldsymbol{\pi}\) 采样一个组件标签 \(\mathbf{z}_n\)。
- 给定 \(z_{nk} = 1\),从伯努利分布 \(\prod_{d=1}^D \mu_{kd}^{x_{nd}} (1 - \mu_{kd})^{1 - x_{nd}}\) 中采样出观测向量 \(\mathbf{x}_n\)。
2. 完全数据的对数似然
- 完全数据是指同时包含了观测 \(\mathbf{x}_n\) 和隐变量 \(\mathbf{z}_n\) 的数据。其联合概率分布为:
\[ p(\mathbf{x}_n, \mathbf{z}_n | \boldsymbol{\pi}, \boldsymbol{\mu}) = \prod_{k=1}^K \left[ \pi_k \prod_{d=1}^D \mu_{kd}^{x_{nd}} (1 - \mu_{kd})^{1 - x_{nd}} \right]^{z_{nk}}. \]
- 对于所有 \(N\) 个独立同分布的数据点,完全数据的对数似然为:
\[ \log p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\pi}, \boldsymbol{\mu}) = \sum_{n=1}^N \sum_{k=1}^K z_{nk} \left[ \log \pi_k + \sum_{d=1}^D \left( x_{nd} \log \mu_{kd} + (1 - x_{nd}) \log (1 - \mu_{kd}) \right) \right]. \]
3. E步:计算隐变量的后验期望(责任值)
- 在EM算法的E步,我们需要在给定当前参数估计 \(\boldsymbol{\pi}^{\text{old}}, \boldsymbol{\mu}^{\text{old}}\) 和观测数据 \(\mathbf{X}\) 的条件下,计算隐变量 \(z_{nk}\) 的后验期望,这个期望通常称为“责任值”(responsibility),记作 \(\gamma_{nk}\)。
- 根据贝叶斯定理,责任值的计算为:
\[ \gamma_{nk} = \mathbb{E}[z_{nk} | \mathbf{x}_n, \boldsymbol{\pi}^{\text{old}}, \boldsymbol{\mu}^{\text{old}}] = p(z_{nk}=1 | \mathbf{x}_n, \boldsymbol{\pi}^{\text{old}}, \boldsymbol{\mu}^{\text{old}}). \]
- 具体地,先计算第 \(k\) 个组件对数据点 \(n\) 的“非归一化后验”:
\[ \tilde{\gamma}_{nk} = \pi_k^{\text{old}} \prod_{d=1}^D (\mu_{kd}^{\text{old}})^{x_{nd}} (1 - \mu_{kd}^{\text{old}})^{1 - x_{nd}}. \]
- 然后进行归一化,得到责任值:
\[ \gamma_{nk} = \frac{\tilde{\gamma}_{nk}}{\sum_{j=1}^K \tilde{\gamma}_{nj}}. \]
- 直观上,\(\gamma_{nk}\) 衡量了数据点 \(n\) 由第 \(k\) 个组件生成的可能性大小,它是一个介于0和1之间的数,且对每个 \(n\) 满足 \(\sum_{k=1}^K \gamma_{nk} = 1\)。
4. M步:最大化期望似然,更新模型参数
- 在M步,我们需要最大化关于参数 \(\boldsymbol{\pi}, \boldsymbol{\mu}\) 的“期望完全数据对数似然”(即Q函数):
\[ Q(\boldsymbol{\pi}, \boldsymbol{\mu}; \boldsymbol{\pi}^{\text{old}}, \boldsymbol{\mu}^{\text{old}}) = \mathbb{E}_{\mathbf{Z} | \mathbf{X}, \boldsymbol{\pi}^{\text{old}}, \boldsymbol{\mu}^{\text{old}}} \left[ \log p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\pi}, \boldsymbol{\mu}) \right]. \]
- 将E步计算出的责任值 \(\gamma_{nk}\) 代入期望,得到:
\[ Q = \sum_{n=1}^N \sum_{k=1}^K \gamma_{nk} \left[ \log \pi_k + \sum_{d=1}^D \left( x_{nd} \log \mu_{kd} + (1 - x_{nd}) \log (1 - \mu_{kd}) \right) \right]. \]
-
现在,我们分别对 \(\boldsymbol{\pi}\) 和 \(\boldsymbol{\mu}\) 最大化 \(Q\) 函数。
(a) 更新混合权重 \(\pi_k\)
- 在约束 \(\sum_{k=1}^K \pi_k = 1\) 下最大化 \(\sum_{n=1}^N \sum_{k=1}^K \gamma_{nk} \log \pi_k\)。利用拉格朗日乘子法,构造拉格朗日函数:
\[ \mathcal{L}_\pi = \sum_{n=1}^N \sum_{k=1}^K \gamma_{nk} \log \pi_k + \lambda \left( \sum_{k=1}^K \pi_k - 1 \right). \]
- 对 \(\pi_k\) 求偏导并令其为零:
\[ \frac{\partial \mathcal{L}_\pi}{\partial \pi_k} = \frac{1}{\pi_k} \sum_{n=1}^N \gamma_{nk} + \lambda = 0 \quad \Rightarrow \quad \sum_{n=1}^N \gamma_{nk} = -\lambda \pi_k. \]
- 对所有 \(k\) 求和,利用 \(\sum_{k=1}^K \pi_k = 1\),可得 \(\sum_{k=1}^K \sum_{n=1}^N \gamma_{nk} = -\lambda\),而左边 \(\sum_{k=1}^K \sum_{n=1}^N \gamma_{nk} = N\)。因此 \(-\lambda = N\)。
- 代入前式,得到更新公式:
\[ \pi_k^{\text{new}} = \frac{1}{N} \sum_{n=1}^N \gamma_{nk}. \]
直观解释:新的混合权重 $ \pi_k $ 等于所有数据点对第 $ k $ 个组件的平均责任。
(b) 更新伯努利参数 \(\mu_{kd}\)
- 在 \(Q\) 函数中,与 \(\mu_{kd}\) 相关的项为:
\[ \sum_{n=1}^N \gamma_{nk} \left[ x_{nd} \log \mu_{kd} + (1 - x_{nd}) \log (1 - \mu_{kd}) \right]. \]
- 对 \(\mu_{kd}\) 求偏导(注意 \(\mu_{kd}\) 只出现在该项中)并令其为零:
\[ \frac{\partial}{\partial \mu_{kd}} = \sum_{n=1}^N \gamma_{nk} \left[ \frac{x_{nd}}{\mu_{kd}} - \frac{1 - x_{nd}}{1 - \mu_{kd}} \right] = 0. \]
- 整理得:
\[ \sum_{n=1}^N \gamma_{nk} x_{nd} \cdot (1 - \mu_{kd}) = \sum_{n=1}^N \gamma_{nk} (1 - x_{nd}) \cdot \mu_{kd}. \]
- 解出 \(\mu_{kd}\):
\[ \mu_{kd}^{\text{new}} = \frac{\sum_{n=1}^N \gamma_{nk} x_{nd}}{\sum_{n=1}^N \gamma_{nk}}. \]
直观解释:新的伯努利参数 $ \mu_{kd} $ 等于所有数据点在第 $ d $ 维特征上取值为1的加权平均,权重为数据点对第 $ k $ 个组件的责任 $ \gamma_{nk} $。
5. 算法流程总结
伯努利混合模型的EM算法可归纳为以下迭代步骤:
- 初始化:随机或启发式地初始化参数 \(\boldsymbol{\pi}^{(0)}, \boldsymbol{\mu}^{(0)}\),设置收敛阈值 \(\epsilon\) 和最大迭代次数。
- 循环直到收敛(例如,参数变化小于 \(\epsilon\) 或达到最大迭代次数):
- E步:用当前参数 \(\boldsymbol{\pi}^{(t)}, \boldsymbol{\mu}^{(t)}\) 计算责任值 \(\gamma_{nk}^{(t)}\):
\[ \gamma_{nk}^{(t)} = \frac{\pi_k^{(t)} \prod_{d=1}^D (\mu_{kd}^{(t)})^{x_{nd}} (1 - \mu_{kd}^{(t)})^{1 - x_{nd}}}{\sum_{j=1}^K \pi_j^{(t)} \prod_{d=1}^D (\mu_{jd}^{(t)})^{x_{nd}} (1 - \mu_{jd}^{(t)})^{1 - x_{nd}}}. \]
- M步:用计算出的 \(\gamma_{nk}^{(t)}\) 更新参数:
\[ \pi_k^{(t+1)} = \frac{1}{N} \sum_{n=1}^N \gamma_{nk}^{(t)}, \]
\[ \mu_{kd}^{(t+1)} = \frac{\sum_{n=1}^N \gamma_{nk}^{(t)} x_{nd}}{\sum_{n=1}^N \gamma_{nk}^{(t)}}. \]
- 输出:最终的参数估计 \(\boldsymbol{\pi}^*, \boldsymbol{\mu}^*\),以及每个数据点的软聚类标签(即责任值 \(\gamma_{nk}\))。
关键要点:
- EM算法通过引入隐变量的后验分布(E步),将复杂的最大似然估计问题转化为一系列简单的、可解析求解的更新(M步)。
- 在伯努利混合模型中,E步计算责任值本质上是进行一个“软分配”,M步更新参数则类似于在“加权数据”上计算伯努利分布的极大似然估计。
- 该模型可广泛应用于二值数据的聚类(如文档的主题建模中的词项出现与否、黑白图像建模等)。