期望最大化（EM）算法在混合伯努利模型（Mixture of Bernoulli）中的参数估计过程
题目描述
假设我们有一组二值数据（每个特征取值为0或1），例如多个文档的词袋表示（每个词出现与否）、多个用户的二进制行为记录等。我们假设这些数据是由一个混合伯努利模型生成的：
- 存在 \(K\) 个伯努利分布组件（即“主题”或“类别”）。
- 每个伯努利组件 \(k\) 由一个 \(D\) 维参数向量 \(\boldsymbol{\mu}_k = (\mu_{k1}, \dots, \mu_{kD})\) 刻画，其中 \(\mu_{kd} \in (0,1)\) 表示在第 \(k\) 个组件中，第 \(d\) 个特征取1的概率。
- 每个数据点 \(\boldsymbol{x}_n = (x_{n1}, \dots, x_{nD}) \in \{0,1\}^D\) 的生成过程为：
- 先选择一个隐变量（组件指示） \(z_n \in \{1,\dots,K\}\)，其概率为混合权重 \(\pi_k = P(z_n = k)\)，满足 \(\sum_{k=1}^K \pi_k = 1\)。
- 给定 \(z_n = k\)，\(\boldsymbol{x}_n\) 的每个维度独立地从伯努利分布生成：
\[ P(x_{nd} = 1 \mid z_n = k) = \mu_{kd}, \quad P(x_{nd} = 0 \mid z_n = k) = 1 - \mu_{kd}. \]
即 $p(\boldsymbol{x}_n \mid z_n = k) = \prod_{d=1}^D \mu_{kd}^{x_{nd}} (1 - \mu_{kd})^{1 - x_{nd}}$。
目标：在仅观测到数据集 \(\mathcal{D} = \{\boldsymbol{x}_1, \dots, \boldsymbol{x}_N\}\) 的情况下，估计模型参数 \(\boldsymbol{\theta} = \{\pi_k, \boldsymbol{\mu}_k\}_{k=1}^K\)。由于存在隐变量 \(z_n\)，直接极大化似然函数困难，因此使用期望最大化（EM）算法迭代求解。
EM算法的整体思路
EM算法是一种迭代优化方法，用于含有隐变量的概率模型的极大似然估计。每次迭代分为两步：
- E步（Expectation）：基于当前参数 \(\boldsymbol{\theta}^{old}\)，计算隐变量 \(z_n\) 的后验分布 \(P(z_n = k \mid \boldsymbol{x}_n, \boldsymbol{\theta}^{old})\)，进而计算完全数据的对数似然函数关于隐变量后验的期望，称为Q函数。
- M步（Maximization）：最大化Q函数，得到新的参数估计 \(\boldsymbol{\theta}^{new}\)。
重复E步和M步直至参数收敛。
第一步：写出完全数据的对数似然函数
完全数据包括观测 \(\boldsymbol{x}_n\) 和隐变量 \(z_n\)。
- 隐变量 \(z_n\) 的分布： \(P(z_n = k) = \pi_k\)。
- 给定 \(z_n = k\)，观测 \(\boldsymbol{x}_n\) 的条件概率（伯努利乘积）：
\[ P(\boldsymbol{x}_n \mid z_n = k, \boldsymbol{\mu}_k) = \prod_{d=1}^D \mu_{kd}^{x_{nd}} (1 - \mu_{kd})^{1 - x_{nd}}. \]
因此，完全数据的似然函数为：
\[P(\boldsymbol{x}_n, z_n = k \mid \boldsymbol{\theta}) = \pi_k \prod_{d=1}^D \mu_{kd}^{x_{nd}} (1 - \mu_{kd})^{1 - x_{nd}}. \]
完全数据的对数似然（所有数据点独立同分布）为：
\[\log P(\{\boldsymbol{x}_n, z_n\} \mid \boldsymbol{\theta}) = \sum_{n=1}^N \left[ \log \pi_{z_n} + \sum_{d=1}^D \big( x_{nd} \log \mu_{z_n, d} + (1 - x_{nd}) \log (1 - \mu_{z_n, d}) \big) \right]. \]
第二步：E步——计算隐变量的后验概率（责任值）
给定当前参数 \(\boldsymbol{\theta}^{old} = \{\pi_k^{old}, \boldsymbol{\mu}_k^{old}\}\)，计算每个数据点 \(\boldsymbol{x}_n\) 属于每个组件 \(k\) 的后验概率，也称为责任值 \(\gamma_{nk}\)：
\[\gamma_{nk} = P(z_n = k \mid \boldsymbol{x}_n, \boldsymbol{\theta}^{old}) = \frac{\pi_k^{old} \, P(\boldsymbol{x}_n \mid z_n = k, \boldsymbol{\mu}_k^{old})}{\sum_{j=1}^K \pi_j^{old} \, P(\boldsymbol{x}_n \mid z_n = j, \boldsymbol{\mu}_j^{old})}. \]
其中，
\[P(\boldsymbol{x}_n \mid z_n = k, \boldsymbol{\mu}_k^{old}) = \prod_{d=1}^D (\mu_{kd}^{old})^{x_{nd}} (1 - \mu_{kd}^{old})^{1 - x_{nd}}. \]
计算技巧：实际计算时，为避免数值下溢，通常对分子和分母同时取对数，在log空间计算，再通过指数化和归一化得到 \(\gamma_{nk}\)。
第三步：构造Q函数
Q函数是完全数据对数似然关于隐变量后验分布的期望：
\[Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{old}) = \mathbb{E}_{P(\{z_n\} \mid \{\boldsymbol{x}_n\}, \boldsymbol{\theta}^{old})} \left[ \log P(\{\boldsymbol{x}_n, z_n\} \mid \boldsymbol{\theta}) \right]. \]
由于各数据点独立，且隐变量后验可分解为 \(\prod_{n=1}^N P(z_n \mid \boldsymbol{x}_n, \boldsymbol{\theta}^{old})\)，我们有：
\[Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{old}) = \sum_{n=1}^N \sum_{k=1}^K \gamma_{nk} \left[ \log \pi_k + \sum_{d=1}^D \big( x_{nd} \log \mu_{kd} + (1 - x_{nd}) \log (1 - \mu_{kd}) \big) \right]. \]
第四步：M步——最大化Q函数以更新参数
我们需要在约束 \(\sum_{k=1}^K \pi_k = 1\) 和 \(0 < \mu_{kd} < 1\) 下，最大化 \(Q\)。
1. 更新混合权重 \(\pi_k\)
引入拉格朗日乘子 \(\lambda\) 处理约束 \(\sum_k \pi_k = 1\)，构造拉格朗日函数：
\[\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 \sum_{n=1}^N \gamma_{nk} = -\lambda \sum_{k=1}^K \pi_k \quad \Rightarrow \quad N = -\lambda \quad \text{（因为 } \sum_{k=1}^K \gamma_{nk} = 1 \text{）}. \]
因此，
\[\pi_k^{new} = \frac{1}{N} \sum_{n=1}^N \gamma_{nk}. \]
直观解释：新的 \(\pi_k\) 是所有数据点属于组件 \(k\) 的后验概率的平均。
2. 更新伯努利参数 \(\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]. \]
令导数为零：
\[\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} (1 - \mu_{kd}) = \sum_{n=1}^N \gamma_{nk} (1 - x_{nd}) \mu_{kd}. \]
进一步：
\[\sum_{n=1}^N \gamma_{nk} x_{nd} - \mu_{kd} \sum_{n=1}^N \gamma_{nk} x_{nd} = \mu_{kd} \sum_{n=1}^N \gamma_{nk} - \mu_{kd} \sum_{n=1}^N \gamma_{nk} x_{nd}. \]
化简后：
\[\sum_{n=1}^N \gamma_{nk} x_{nd} = \mu_{kd} \sum_{n=1}^N \gamma_{nk}. \]
因此，
\[\mu_{kd}^{new} = \frac{\sum_{n=1}^N \gamma_{nk} x_{nd}}{\sum_{n=1}^N \gamma_{nk}}. \]
直观解释：新的 \(\mu_{kd}\) 是所有数据点在第 \(d\) 维上取1的加权平均，权重为数据点属于组件 \(k\) 的后验概率 \(\gamma_{nk}\)。
第五步：迭代执行E步和M步
- 初始化参数：随机或启发式设置 \(\pi_k^{(0)}\) 和 \(\mu_{kd}^{(0)}\)，通常 \(\pi_k^{(0)} = 1/K\)，\(\mu_{kd}^{(0)}\) 在 \((0,1)\) 内随机生成。
- 重复以下步骤直至收敛（参数变化小于阈值或似然函数不再显著增加）：
- E步：用当前参数 \(\boldsymbol{\theta}^{(t)}\) 计算责任值 \(\gamma_{nk}^{(t+1)}\)。
- M步：用 \(\gamma_{nk}^{(t+1)}\) 按上述公式计算 \(\pi_k^{(t+1)}\) 和 \(\mu_{kd}^{(t+1)}\)。
收敛性：EM算法保证每次迭代后不完全数据的对数似然 \(P(\{\boldsymbol{x}_n\} \mid \boldsymbol{\theta})\) 不会减少，最终收敛到局部极大值。
总结
- 混合伯努利模型适用于二值数据的软聚类（每个点以概率属于各个组件）。
- EM算法的核心：通过引入隐变量 \(z_n\)，将极大似然估计转化为迭代的E步（计算后验责任）和M步（更新参数）。
- E步：计算 \(\gamma_{nk} = P(z_n = k \mid \boldsymbol{x}_n, \boldsymbol{\theta}^{old})\)。
- M步：更新 \(\pi_k = \frac{1}{N} \sum_n \gamma_{nk}\)，\(\mu_{kd} = \frac{\sum_n \gamma_{nk} x_{nd}}{\sum_n \gamma_{nk}}\)。
这个算法在文本聚类（每个文档表示为词频二值向量）、二值图像建模等场景有直接应用。