EM算法 Python 3.12 实现：硬币实验单次迭代收敛速度实测（附完整代码）

EM算法Python实现：硬币实验单次迭代收敛速度深度解析

1. EM算法核心思想与硬币实验场景

EM算法作为机器学习经典方法，其核心在于通过E步（期望计算）和M步（最大化）的交替迭代，解决含隐变量的概率模型参数估计问题。硬币实验作为经典示例，完美展示了EM算法的运作机制：

• 实验设定：假设有两枚质地不同的硬币A和B，每次实验随机选择一枚进行多次抛掷
• 观测数据：记录每次抛掷结果（正面1/反面0），但不记录使用的是哪枚硬币
• 核心挑战：在不知道每次实验所用硬币的情况下，估计两枚硬币各自的正面概率

import numpy as np from scipy import stats def simulate_coin_toss(pA=0.4, pB=0.6, n_experiments=10, n_tosses=5): """生成模拟硬币实验数据""" choices = np.random.choice(['A','B'], size=n_experiments) observations = [] for coin in choices: p = pA if coin == 'A' else pB obs = np.random.binomial(1, p, size=n_tosses) observations.append(obs.tolist()) return observations

2. 单次迭代的数学原理与实现

2.1 E步骤：隐变量概率估计

在E步骤中，我们基于当前参数θ计算隐变量（硬币选择）的后验概率。对于每次实验observation：

1. 计算当前参数下各硬币产生该结果的概率
2. 通过贝叶斯定理得到权重分配

def e_step(observation, theta_A, theta_B): len_obs = len(observation) num_heads = sum(observation) num_tails = len_obs - num_heads # 计算两枚硬币产生该结果的概率 prob_A = stats.binom.pmf(num_heads, len_obs, theta_A) prob_B = stats.binom.pmf(num_heads, len_obs, theta_B) # 归一化得到权重 weight_A = prob_A / (prob_A + prob_B) weight_B = 1 - weight_A return weight_A, weight_B

2.2 M步骤：参数最大化

在M步骤中，我们基于E步得到的权重重新估计参数：

1. 计算各硬币的期望正反面次数
2. 通过极大似然估计更新参数

def m_step(observations, theta_A, theta_B): counts = {'A': {'H': 0, 'T': 0}, 'B': {'H': 0, 'T': 0}} for obs in observations: weight_A, weight_B = e_step(obs, theta_A, theta_B) num_heads = sum(obs) num_tails = len(obs) - num_heads # 更新期望计数 counts['A']['H'] += weight_A * num_heads counts['A']['T'] += weight_A * num_tails counts['B']['H'] += weight_B * num_heads counts['B']['T'] += weight_B * num_tails # 计算新参数 new_theta_A = counts['A']['H'] / (counts['A']['H'] + counts['A']['T']) new_theta_B = counts['B']['H'] / (counts['B']['H'] + counts['B']['T']) return new_theta_A, new_theta_B

3. 完整EM算法实现与收敛分析

3.1 完整迭代流程

将E步和M步结合，实现完整的EM算法：

def em_algorithm(observations, initial_theta, tol=1e-6, max_iter=100): theta_A, theta_B = initial_theta history = [initial_theta] for i in range(max_iter): # M步 new_theta_A, new_theta_B = m_step(observations, theta_A, theta_B) # 检查收敛 delta = abs(new_theta_A - theta_A) + abs(new_theta_B - theta_B) if delta < tol: break theta_A, theta_B = new_theta_A, new_theta_B history.append((theta_A, theta_B)) return (theta_A, theta_B), history

3.2 收敛速度实测

我们通过实验分析不同初始值对收敛速度的影响：

注意：实验结果基于模拟数据，真实值θA=0.4，θB=0.6。初始值接近真实值时收敛更快。

4. 可视化分析与性能优化

4.1 收敛过程可视化

import matplotlib.pyplot as plt def plot_convergence(history): plt.figure(figsize=(10, 6)) theta_A = [x[0] for x in history] theta_B = [x[1] for x in history] plt.plot(theta_A, label='θA', marker='o') plt.plot(theta_B, label='θB', marker='s') plt.axhline(0.4, color='red', linestyle='--', alpha=0.3) plt.axhline(0.6, color='blue', linestyle='--', alpha=0.3) plt.xlabel('Iteration') plt.ylabel('Parameter Value') plt.title('EM Algorithm Convergence') plt.legend() plt.grid(True) plt.show()

4.2 数值稳定性优化

实际实现中需注意数值稳定性问题：

def stable_e_step(observation, theta_A, theta_B, epsilon=1e-10): len_obs = len(observation) num_heads = sum(observation) # 添加极小值避免零概率 prob_A = stats.binom.pmf(num_heads, len_obs, theta_A) + epsilon prob_B = stats.binom.pmf(num_heads, len_obs, theta_B) + epsilon # 对数空间计算提高数值稳定性 log_prob_A = np.log(prob_A) log_prob_B = np.log(prob_B) max_log = max(log_prob_A, log_prob_B) weight_A = np.exp(log_prob_A - max_log) weight_B = np.exp(log_prob_B - max_log) # 归一化 total = weight_A + weight_B return weight_A / total, weight_B / total

5. 工程实践中的关键考量

1. 初始值选择：实践中建议：

   • 运行多次EM算法，选择不同随机初始值
   • 选择似然函数值最大的结果作为最终解

2. 停止准则：除参数变化外，还可监测：

   • 对数似然函数的变化量
   • 最大迭代次数的合理设置

3. 高维扩展：当处理更复杂模型时：

   • 考虑使用加速EM算法变种
   • 并行化E步骤计算

def log_likelihood(observations, theta_A, theta_B): total = 0.0 for obs in observations: num_heads = sum(obs) len_obs = len(obs) # 混合概率 prob_A = stats.binom.pmf(num_heads, len_obs, theta_A) prob_B = stats.binom.pmf(num_heads, len_obs, theta_B) total += np.log(0.5 * prob_A + 0.5 * prob_B + 1e-10) return total

硬币实验虽然简单，但完整展现了EM算法的核心思想。在实际项目中遇到更复杂的隐变量模型时，这个实现框架仍具有指导意义。理解这个基础案例后，可以更容易地将其扩展到高斯混合模型、隐马尔可夫模型等更复杂的场景。