期望最大化（EM）算法在混合泊松模型（Mixture of Poisson Models）中的参数估计过程
我将为你详细讲解EM算法在混合泊松模型中的参数估计过程，包括背景知识、问题定义、数学推导和计算步骤。
一、问题背景介绍
混合泊松模型是计数数据（如网站访问量、客服中心电话数、交通事故次数等）分析中常用的模型。假设我们有观测到的计数数据，但不知道每个数据点来自哪个具体的泊松分布。通过EM算法，我们可以同时估计：
- 每个泊松分布的参数（强度参数λ）
- 每个数据点属于每个分布的"软"分配（概率权重）
二、混合泊松模型的形式化定义
假设我们有K个泊松分布混合，每个分布的强度参数为λ_k（k=1,2,...,K），混合权重为π_k（满足∑π_k=1且π_k≥0）。
对于第i个观测数据x_i，其概率为：
P(x_i|Θ) = ∑_{k=1}^K π_k * P(x_i|λ_k)
其中Θ = {π_k, λ_k} 是所有参数，P(x_i|λ_k)是泊松分布的概率质量函数：
P(x_i|λ_k) = (λ_k^{x_i} * e^{-λ_k}) / x_i!
三、EM算法的E步骤：计算隐变量的期望
-
隐变量定义：
定义隐变量z_{ik}，当观测i属于第k个泊松分布时z_{ik}=1，否则为0。实际上z_{ik}是未知的。 -
计算后验概率（责任γ_{ik}）：
在给定当前参数估计的条件下，计算观测i来自第k个分量的后验概率：
γ_{ik} = P(z_{ik}=1|x_i, Θ^{(t)})
= [π_k^{(t)} * P(x_i|λ_k^{(t)})] / ∑_{j=1}^K [π_j^{(t)} * P(x_i|λ_j^{(t)})]
这里：
- π_k^{(t)}是第t轮迭代时第k个分量的混合权重
- λ_k^{(t)}是第t轮迭代时第k个分量的泊松强度参数
- P(x_i|λ_k^{(t)}) = (λ_k^{(t)})^{x_i} * e^{-λ_k^{(t)}} / x_i!
- 几何解释：γ_{ik}是"软"分配概率，表示第i个数据点属于第k个分量的置信度。所有K个分量的γ_{ik}之和为1。
四、EM算法的M步骤：最大化期望
-
目标函数：完全数据的对数似然函数在给定当前参数下的期望
Q(Θ|Θ^{(t)}) = E_{Z|X,Θ^{(t)}}[log P(X,Z|Θ)]
= ∑{i=1}^N ∑{k=1}^K γ_{ik} * [log π_k + log P(x_i|λ_k)] -
更新混合权重π_k：
在约束∑π_k=1下最大化Q函数：
π_k^{(t+1)} = (1/N) * ∑{i=1}^N γ{ik}
直观解释：第k个分量的新权重是"软"分配给它的所有数据点比例的平均。 -
更新泊松强度参数λ_k：
对每个分量k，令∂Q/∂λ_k = 0：
∑{i=1}^N γ{ik} * [∂/∂λ_k log P(x_i|λ_k)] = 0
其中log P(x_i|λ_k) = x_i * log λ_k - λ_k - log(x_i!)
求偏导：
∂/∂λ_k [x_i * log λ_k - λ_k] = x_i/λ_k - 1
设偏导为0：
∑{i=1}^N γ{ik} * (x_i/λ_k - 1) = 0
=> λ_k = [∑{i=1}^N γ{ik} * x_i] / ∑{i=1}^N γ{ik}
因此：
λ_k^{(t+1)} = [∑{i=1}^N γ{ik} * x_i] / ∑{i=1}^N γ{ik}
直观解释：第k个分量的新强度参数是分配给它的所有数据点的加权平均计数，权重是责任γ_{ik}。
五、完整的EM算法迭代步骤
步骤1：参数初始化
- 随机初始化混合权重π_k^{(0)}（需满足和为1）
- 随机初始化强度参数λ_k^{(0)}，通常设为均匀分布或基于数据范围
- 设置收敛阈值ε和最大迭代次数T
步骤2：E步骤（第t轮）
对于每个数据点i=1,...,N和每个分量k=1,...,K：
- 计算每个分量的似然：P(x_i|λ_k^{(t)}) = (λ_k^{(t)})^{x_i} * e^{-λ_k^{(t)}} / x_i!
- 计算未归一化的权重：w_{ik} = π_k^{(t)} * P(x_i|λ_k^{(t)})
- 归一化得到责任：γ_{ik} = w_{ik} / ∑{j=1}^K w{ij}
步骤3：M步骤（第t轮）
- 更新混合权重：π_k^{(t+1)} = (1/N) * ∑{i=1}^N γ{ik}
- 更新强度参数：λ_k^{(t+1)} = [∑{i=1}^N γ{ik} * x_i] / ∑{i=1}^N γ{ik}
步骤4：检查收敛
- 计算对数似然的变化：ΔL = |L(Θ^{(t+1)}) - L(Θ^{(t)})|
其中L(Θ) = ∑{i=1}^N log[∑{k=1}^K π_k * P(x_i|λ_k)] - 如果ΔL < ε或达到最大迭代次数T，则停止；否则返回步骤2
六、算法特性与注意事项
-
泊松分布的特性：λ必须为正，这由更新公式自动保证（分子为加权计数和，分母为正）
-
初始化的敏感性：
- λ_k的初始值最好分散在数据范围内
- 多次随机初始化可避免局部最优
- 可用K-means结果作为初始化（但K-means用于计数数据需谨慎）
-
零计数的处理：
当x_i=0时，P(x_i|λ_k) = e^{-λ_k}，计算不会出问题
但需注意数值稳定性，特别当λ_k很小时 -
过拟合问题：
- 可通过交叉验证选择K值
- 或使用信息准则（AIC/BIC）进行模型选择
- 边界情况处理：
- 如果∑γ_{ik}接近0，对应分量可移除
- 数值计算时使用对数似然避免下溢
七、实际应用示例
假设我们要对某客服中心每小时接到的电话数建模，数据明显呈多峰分布（例如：工作时间、午餐时间、夜间等模式不同）。通过混合泊松模型：
- λ_1≈5（夜间低强度）
- λ_2≈20（午餐时间中等强度）
- λ_3≈50（工作时间高强度）
EM算法可自动从数据中学习这些参数以及每个时间段属于哪种模式的比例。
这个算法框架可扩展为更复杂的模型，如加入协变量的混合泊松回归，但核心的EM迭代思路保持不变。