EM算法到底是什么东东

EM（Expectation-Maximization期望最大化）算法是机器学习中非常重要的一类算法，广泛应用于聚类、缺失数据建模、隐变量模型学习等场景，比如高斯混合模型（GMM）就是经典应用。

🐤 第一步：直观理解

EM算法的核心是：

我不知道这个数据是哪一类（隐变量），就先猜；然后根据可见的情况，慢慢猜的更准。
EM算法就是一个“猜→修正→再猜”的循环。

例子1：

• 给你一篇文章让你读
• 可观测数据：文档中的词语。
• 隐变量：文档的主题分布。
• 本质：主题是潜在的，决定了词语的出现概率。

例子2：
假设有两个数据分布（两类），然后随机从这两个分布里抽出一些样本交给你，你不知道给你的样本点属于哪一类（隐含的类别），以及这两个数据分布的统计特性（均值，方差）

EM算法的做法是：

1. 随便猜一下每个点属于哪个类别（初始猜测）
2. 计算：在当前参数下，每个点属于各个类别的“概率”（这是E步）
3. 用这些概率来“反推”出最合理的类别参数（比如均值、方差）（这是M步）
4. 重复步骤2-3，直到参数不怎么变为止。

---

✍️ 第二步：数学公式

你有一堆数据点 ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$，你相信这些数据来自 $K$ 个不同的高斯分布：

• 每个分布 $k$ 有自己的参数：均值 ${\mu }_k$、方差 ${\sigma }_k^2$、权重 ${\pi }_k$（概率总和为1）
• 但你不知道哪个点来自哪个分布（这是隐变量）

---

E步（Expectation），即“先猜”

初始化：随机初始化均值 ${\mu }_k$、方差 ${\sigma }_k^2$ 和权重 ${\pi }_k$

计算每个样本属于每个高斯分布的“后验概率”：

$${\gamma }_{ik}=\frac{{\pi }_k\cdot \mathcal{N}(x_i|{\mu }_k,{\sigma }_k^2)}{{\sum }_{j=1}^K{\pi }_j\cdot \mathcal{N}(x_i|{\mu }_j,{\sigma }_j^2)}$$

这表示：样本 $x_i$ 属于第 $k$ 个高斯分布的概率。

---

M步（Maximization），即“反推参数”

根据这些概率 ${\gamma }_{ik}$ 来重新估计参数：

$${\mu }_k=\frac{{\sum }_i{\gamma }_{ik}{x}_i}{{\sum }_i{\gamma }_{ik}},{\sigma }_k^2=\frac{{\sum }_i{\gamma }_{ik}(x_i-{\mu }_k{)}^2}{{\sum }_i{\gamma }_{ik}},{\pi }_k=\frac{1}{n}\sum _i{\gamma }_{ik}$$

---

🧊 第三步 ：一个具体的例子——高斯混合模型（GMM）

什么是GMM？

高斯混合模型（GMM）就是用多个“高斯分布”加权叠加来组合描述一个复杂的数据分布。GMM 的参数（每个高斯的均值、方差、权重）不能直接算出来，但可以用 EM算法 来一步步逼近！

• GMM = 模型框架
• EM = 参数求解方法

GMM分布可视化

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm# 解决中文显示问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题# 定义两个高斯分布的参数
# 每个分布由均值(mu)、标准差(sigma)和权重(weight)组成
mu1, sigma1, weight1 = 5, 1, 0.4  # 分布1: 均值5, 标准差1, 权重40%
mu2, sigma2, weight2 = 15, 2, 0.6  # 分布2: 均值15, 标准差2, 权重60%# 生成X轴范围，覆盖两个分布的3σ范围
x_min = min(mu1 - 3*sigma1, mu2 - 3*sigma2)
x_max = max(mu1 + 3*sigma1, mu2 + 3*sigma2)
x = np.linspace(x_min, x_max, 1000)  # 在合理范围内生成1000个点# 计算单个分布的概率密度函数(PDF)
pdf1 = weight1 * norm.pdf(x, mu1, sigma1)  # 第一个高斯分布的加权PDF
pdf2 = weight2 * norm.pdf(x, mu2, sigma2)  # 第二个高斯分布的加权PDF# 计算混合后的整体分布(GMM的概率密度)
pdf_total = pdf1 + pdf2  # 高斯混合模型的PDF是两个加权高斯分布的和# 创建图形并设置大小
plt.figure(figsize=(10, 6))# 绘制各个分布
plt.plot(x, pdf1, label=f"高斯分布1 (μ={mu1}, σ={sigma1}, 权重={weight1})",linestyle='--', color='blue')
plt.plot(x, pdf2, label=f"高斯分布2 (μ={mu2}, σ={sigma2}, 权重={weight2})",linestyle='--', color='green')
plt.plot(x, pdf_total, label="混合分布 GMM", linestyle='--',color='red', linewidth=1.5)# 添加图形标题和标签
plt.title("高斯混合模型(GMM)示意图", fontsize=14)
plt.xlabel("特征值 (示例:糖分含量)", fontsize=12)
plt.ylabel("概率密度", fontsize=12)# 添加图例和网格
plt.legend(fontsize=10)
plt.grid(True, linestyle='--', alpha=0.6)# 显示图形
plt.tight_layout()  # 自动调整子图参数，使之填充整个图像区域
plt.show()

问题背景

假设我们有一堆学生身高数据，比如160cm、155cm、175cm等等，但我们不知道每个学生是小学生还是中学生。我们猜测这些身高来自两个群体：

• 小学生：身高服从一个正态分布（高斯分布），有自己的均值和标准差。
• 中学生：身高服从另一个正态分布，也有自己的均值和标准差。

此外，每个群体在总数据中占一定比例。我们的目标是：

1. 弄清楚每个学生属于小学生还是中学生的概率。
2. 估计两个群体的参数：比例（$\pi$）、均值（$\mu$）、标准差（$\sigma$）。

因为我们不知道真实的类别和参数，所以要用EM算法通过迭代来解决这个问题。

---

EM算法是什么？

EM算法（Expectation-Maximization）是一种用来处理“隐变量”问题的工具。这里，隐变量就是“每个学生属于哪个群体”，我们看不到它，但可以通过数据推测。EM算法分为两步：

• E步（期望）：根据当前猜测的参数，算出每个学生属于小学生或中学生的概率。
• M步（最大化）：用这些概率更新参数，让模型更好地拟合数据。

这两步不断重复，直到参数稳定。

---

具体案例：一步步拆解

1. 初始化：随便猜参数

我们先随便猜一下两个群体的参数，作为起点：

• 小学生：
  • 比例（${\pi }_1$）：50%（0.5）
  • 均值（${\mu }_1$）：150cm
  • 标准差（${\sigma }_1$）：5cm
• 中学生：
  • 比例（${\pi }_2$）：50%（0.5）
  • 均值（${\mu }_2$）：170cm
  • 标准差（${\sigma }_2$）：6cm

这些是初始猜测，不一定准确，但EM算法会帮我们调整。

---

2. E步：算概率（责任值）

现在拿一个学生，身高是160cm。我们要算他属于小学生还是中学生的概率。

（1）用正态分布公式算“可能性”

每个群体都有一个正态分布曲线：

• 小学生：均值150cm，标准差5cm。
• 中学生：均值170cm，标准差6cm。

正态分布的公式是：
$P(X)=\frac{1}{\sqrt{2\pi }\cdot \sigma }\exp \left(-\frac{(X-\mu {)}^2}{2{\sigma }^2}\right)$

• 小学生：
  $P(\text{小学生}|160)=\frac{1}{\sqrt{2\pi }\cdot 5}\exp \left(-\frac{(160-150{)}^2}{2\cdot 5^2}\right)$
  计算指数部分：(160 - 150)² = 100，2 × 5² = 50，-100 / 50 = -2，exp(-2) ≈ 0.135。所以结果是一个较小的数。

• 中学生：
  $P(\text{中学生}|160)=\frac{1}{\sqrt{2\pi }\cdot 6}\exp \left(-\frac{(160-170{)}^2}{2\cdot 6^2}\right)$
  计算指数部分：(160 - 170)² = 100，2 × 6² = 72，-100 / 72 ≈ -1.39，exp(-1.39) ≈ 0.25。结果比小学生的稍大。

简单来说：

• 160cm离150cm（小学生均值）较远，所以可能性较低。
• 160cm离170cm（中学生均值）较近，所以可能性较高。

（2）结合比例算后验概率（责任值）

光看可能性还不够，还要考虑每个群体占的总比例（${\pi }_1=0.5$，${\pi }_2=0.5$）。用贝叶斯公式：
$P(\text{小学生}|160)=\frac{{\pi }_1\cdot P(160|\text{小学生})}{{\pi }_1\cdot P(160|\text{小学生})+{\pi }_2\cdot P(160|\text{中学生})}$

假设计算后：

• $P(\text{小学生}|160)\approx 0.3$（30%）
• $P(\text{中学生}|160)\approx 0.7$（70%）

意思是：这个160cm的学生有30%概率是小学生，70%概率是中学生。对所有学生都做类似计算。

---

3. M步：更新参数

现在我们用所有学生的概率来调整参数。假设有3个学生：160cm、155cm、175cm，E步算出的概率如下：

身高	P(小学生)	P(中学生)
160cm	0.3	0.7
155cm	0.6	0.4
175cm	0.1	0.9

（1）更新比例（$\pi$）

• 新${\pi }_1=$所有学生属于小学生的概率平均值：
  ${\pi }_1=\frac{0.3+0.6+0.1}{3}=\frac{1.0}{3}\approx 0.33$
• 新${\pi }_2=1-{\pi }_1\approx 0.67$

（2）更新均值（$\mu$）

• 新${\mu }_1=$身高 × 属于小学生的概率的加权平均：
  ${\mu }_1=\frac{(160\cdot 0.3)+(155\cdot 0.6)+(175\cdot 0.1)}{0.3+0.6+0.1}=\frac{48+93+17.5}{1.0}=158.5\text{cm}$
• 新${\mu }_2=$类似计算：
  ${\mu }_2=\frac{(160\cdot 0.7)+(155\cdot 0.4)+(175\cdot 0.9)}{0.7+0.4+0.9}=\frac{112+62+157.5}{2.0}=165.75\text{cm}$

（3）更新标准差（$\sigma$）

• 新${\sigma }_1=$ 身高偏离新均值${\mu }_1$的加权方差：
  ${\sigma }_1=\sqrt{\frac{0.3\cdot (160-158.5{)}^2+0.6\cdot (155-158.5{)}^2+0.1\cdot (175-158.5{)}^2}{1.0}}$
  计算后可能得到一个新值，比如4.8cm。
• 新${\sigma }_2=$类似计算，得到新值，比如5.5cm。

---

4. 重复迭代

用新参数（${\pi }_1=0.33,{\mu }_1=158.5,{\sigma }_1=4.8,{\pi }_2=0.67,{\mu }_2=165.75,{\sigma }_2=5.5$）再跑一遍E步和M步。每轮迭代后，参数会更接近真实值。重复直到参数几乎不变，比如：

• 小学生：${\pi }_1\approx 0.4,{\mu }_1\approx 148\text{cm},{\sigma }_1\approx 4\text{cm}$
• 中学生：${\pi }_2\approx 0.6,{\mu }_2\approx 172\text{cm},{\sigma }_2\approx 5\text{cm}$

这意味着EM算法成功把混合的身高数据分成了两个群体，并估计了它们的特征。

---

📊 第四步：完整Python代码

import numpy as np
from scipy.stats import norm
from sklearn.cluster import KMeans
import seaborn as sns
import matplotlib.pyplot as pltclass GMM_EM:"""高斯混合模型（GMM）的EM算法核心实现 - 用于学生身高分布分析案例背景：假设数据包含两个学生群体的身高数据：1. 小学生：服从N(μ1, σ1²)2. 中学生：服从N(μ2, σ2²)每个群体在总样本中占有比例π目标：通过EM算法估计这两个群体的分布参数（π, μ, σ）"""def __init__(self, n_components=2, max_iter=100, tol=1e-6, random_state=42):"""模型初始化参数:n_components : int, default=2要区分的学生群体数量（默认2类：小学生/中学生）max_iter : int, default=100EM算法最大迭代次数tol : float, default=1e-6参数变化收敛阈值（当参数变化小于此值时停止迭代）random_state : int, default=42随机种子，保证结果可重复"""self.n_components = n_components  # 学生群体数量self.max_iter = max_iter  # 最大迭代次数self.tol = tol  # 收敛判断阈值self.random_state = random_state  # 随机种子self.pi = None  # 各群体比例（小学生/中学生的样本占比）self.mu = None  # 各群体身高均值（单位：厘米）self.sigma = None  # 各群体身高标准差self.converged = False  # 是否收敛标志self.iterations = 0  # 实际迭代次数def _validate_input(self, data):"""输入验证 - 确保数据适合学生身高分析验证条件：1. 输入必须是1维数组（每个元素代表一个学生的身高）2. 样本量必须大于群体数量（防止无法区分群体）"""if not isinstance(data, np.ndarray) or data.ndim != 1:raise ValueError("输入应为1维数组，表示学生身高测量值")if len(data) < self.n_components:raise ValueError("样本量需大于群体数量才能进行有效分析")def _initialize_parameters(self, data):"""参数初始化 - 使用K-means进行初步群体划分初始化策略：1. 通过K-means将学生按身高初步分为n_components个群体2. 按群体身高均值升序排列（确保小学生群体在前）3. 初始化参数：- π: 各群体样本占比- μ: 各群体身高均值- σ: 各群体身高标准差（至少1cm防止数值问题）"""np.random.seed(self.random_state)kmeans = KMeans(n_clusters=self.n_components,random_state=self.random_state)labels = kmeans.fit_predict(data.reshape(-1, 1))# 按身高均值升序排列群体（保证小学生群体在前）unique_labels = np.unique(labels)means = np.array([data[labels == lbl].mean() for lbl in unique_labels])order = np.argsort(means)# 初始化参数self.pi = np.array([np.mean(labels == lbl) for lbl in unique_labels[order]])self.mu = np.array([data[labels == lbl].mean() for lbl in unique_labels[order]])self.sigma = np.array([data[labels == lbl].std() if np.sum(labels == lbl) > 1 else 1.0for lbl in unique_labels[order]])def _e_step(self, data):"""期望步（E-step）- 计算学生归属各群体的后验概率计算公式：P(群体k|身高) = π_k * N(身高|μ_k, σ_k²) / Σ(π_j * N(身高|μ_j, σ_j²))返回:responsibilities : array, shape (n_samples, n_components)每个学生属于各群体的概率矩阵"""# 计算各群体的概率密度pdf = norm.pdf(data[:, np.newaxis], self.mu, self.sigma)# 计算责任矩阵（未归一化的后验概率）responsibilities = self.pi * pdf# 归一化使各学生概率和为1responsibilities /= responsibilities.sum(axis=1, keepdims=True)return responsibilitiesdef _m_step(self, data, responsibilities):"""最大化步（M-step）- 更新群体参数更新公式：1. π_k = 群体k的责任值总和 / 总样本数2. μ_k = Σ(责任值_ki * 身高_i) / 群体k的责任值总和3. σ_k = sqrt(Σ(责任值_ki * (身高_i - μ_k)^2) / 群体k的责任值总和)"""# 各群体有效样本数N_k = responsibilities.sum(axis=0)# 更新群体比例self.pi = N_k / len(data)# 更新群体均值self.mu = np.dot(responsibilities.T, data) / N_k# 更新群体标准差diff_sq = (data[:, np.newaxis] - self.mu) ** 2self.sigma = np.sqrt(np.sum(responsibilities * diff_sq, axis=0) / N_k)def _has_converged(self, prev_params):"""收敛判断 - 检查参数是否稳定判断标准：新旧参数（π, μ, σ）的变化是否均小于tol阈值"""return all(np.allclose(new, old, atol=self.tol) for new, old inzip([self.pi, self.mu, self.sigma], prev_params))def fit(self, data):"""训练模型 - EM算法主循环执行流程：1. 输入验证2. 参数初始化3. 迭代执行E步和M步4. 检查收敛或达到最大迭代次数5. 最终按身高均值排序群体"""self._validate_input(data)self._initialize_parameters(data)self.converged = Falsefor self.iterations in range(1, self.max_iter + 1):prev_params = [arr.copy() for arr in [self.pi, self.mu, self.sigma]]# E步：计算后验概率responsibilities = self._e_step(data)# M步：更新参数self._m_step(data, responsibilities)# 检查收敛if self._has_converged(prev_params):self.converged = Truebreak# 最终按身高均值排序群体（保证小学生群体在前）order = np.argsort(self.mu)self.mu = self.mu[order]self.pi = self.pi[order]self.sigma = self.sigma[order]return selfdef predict_proba(self, data):"""预测概率 - 返回每个学生属于各群体的概率返回:array, shape (n_samples, n_components)每个元素表示对应学生属于该群体的概率"""return self._e_step(data)def predict(self, data):"""预测类别 - 返回最可能的群体标签返回:array, shape (n_samples,)每个元素为0（小学生）或1（中学生）"""return np.argmax(self.predict_proba(data), axis=1)def get_params(self):"""获取训练后的模型参数返回:dict 包含：- pi : 各群体比例- mu : 各群体平均身高（cm）- sigma : 各群体身高标准差（cm）- converged : 是否收敛- iterations : 实际迭代次数"""return {'pi': self.pi,'mu': self.mu,'sigma': self.sigma,'converged': self.converged,'iterations': self.iterations}class GMM_Visualizer:"""GMM可视化工具类"""def __init__(self, model, data, true_labels=None):self.model = modelself.data = dataself.true_labels = true_labelsself.colors = sns.color_palette("husl", model.n_components)self.group_labels = ['Primary school student', 'Middle school student'] \if model.n_components == 2 else [f'Group {i+1}' for i in range(model.n_components)]def plot_results(self, save_path=None):"""可视化拟合结果"""plt.figure(figsize=(12, 7))x = np.linspace(self.data.min()-15, self.data.max()+15, 1000)# 绘制直方图sns.histplot(self.data, bins=30, kde=False, stat='density',color='gray', alpha=0.3, label='Original Data')# 绘制各成分和混合分布mixture_pdf = np.zeros_like(x)for k in range(self.model.n_components):component_pdf = self.model.pi[k] * norm.pdf(x, self.model.mu[k], self.model.sigma[k])plt.plot(x, component_pdf, color=self.colors[k], lw=2,label=f'{self.group_labels[k]} (π={self.model.pi[k]:.2f}, μ={self.model.mu[k]:.1f}, σ={self.model.sigma[k]:.1f})')mixture_pdf += component_pdf# 绘制混合分布plt.plot(x, mixture_pdf, 'k--', lw=2.5, label='Mixture Distribution')# 绘制真实分布（如果存在）if self.true_labels is not None:for k in range(self.model.n_components):sns.histplot(self.data[self.true_labels == k], bins=15, kde=False, stat='density',color=self.colors[k], alpha=0.3, label=f'True {self.group_labels[k]}')plt.title('GMM Fitting Results', fontsize=14)plt.xlabel('Height (cm)')plt.ylabel('Probability Density')plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', frameon=True)plt.grid(alpha=0.2)if save_path:plt.savefig(save_path, dpi=300, bbox_inches='tight')print(f"图像已保存到: {save_path}")else:plt.show()plt.close()def generate_data(n_primary=100, n_secondary=100,primary_mean=160, primary_std=5,secondary_mean=175, secondary_std=7,random_state=42):"""生成模拟学生身高数据案例背景：生成包含两个学生群体的身高数据集，模拟实际观测数据：- 小学生群体：身高服从正态分布N(μ1, σ1²)- 中学生群体：身高服从正态分布N(μ2, σ2²)参数:n_primary : int, default=100小学生样本数量（默认100人）n_secondary : int, default=100中学生样本数量（默认100人）primary_mean : float, default=160小学生群体平均身高（厘米）primary_std : float, default=5小学生群体身高标准差（厘米）secondary_mean : float, default=175中学生群体平均身高（厘米）secondary_std : float, default=7中学生群体身高标准差（厘米）random_state : int, default=42随机种子，保证数据生成可重复返回:data : ndarray, shape (n_samples,)打乱后的混合身高数据（单位：厘米）labels : ndarray, shape (n_samples,)对应的学生群体标签（0:小学生, 1:中学生）"""# 设置随机种子保证结果可重复np.random.seed(random_state)# 生成小学生身高数据（正态分布）primary = np.random.normal(primary_mean, primary_std, n_primary)# 生成中学生身高数据（正态分布）secondary = np.random.normal(secondary_mean, secondary_std, n_secondary)# 合并数据并创建标签data = np.concatenate([primary, secondary])labels = np.concatenate([np.zeros(n_primary), np.ones(n_secondary)])# 打乱数据以模拟真实观测场景# 保持数据与标签的对应关系idx = np.random.permutation(len(data))return data[idx], labels[idx].astype(int)def print_results(model, data, true_labels=None):"""打印结果对比"""params = model.get_params()pi, mu, sigma = params['pi'], params['mu'], params['sigma']print("\n" + "=" * 60)print("GMM-EM 模型预测结果")print("=" * 60)for k in range(len(pi)):print(f"Group {k+1}: 比例={pi[k]:.4f}, 均值={mu[k]:.2f}cm, 方差={sigma[k]:.2f}cm")if true_labels is not None:print("\n真实参数:")for k in range(len(pi)):mask = true_labels == kprint(f"Group {k+1}: 比例={np.mean(mask):.4f}, "f"均值={data[mask].mean():.2f}cm, 方差={data[mask].std():.2f}cm")accuracy = np.mean(model.predict(data) == true_labels)print(f"\n分类准确率: {accuracy:.2%}")print(f"是否收敛: {params['converged']}, 迭代次数: {params['iterations']}")print("=" * 60)if __name__ == "__main__":# 生成数据data, labels = generate_data()# 训练模型model = GMM_EM(n_components=2)model.fit(data)# 打印结果print_results(model, data, labels)# 可视化visualizer = GMM_Visualizer(model, data, labels)visualizer.plot_results(save_path="./fitted_distribution.png")

📌 第五步：总结重点

概念	含义
隐变量	不知道但是存在的变量（比如样本的真实类别）
E步	计算每个数据点属于哪个分布的“概率”
M步	根据这个概率重新计算每个分布的参数
收敛	参数变化很小，不再更新，算法停止
应用场景	高斯混合聚类、缺失数据估计、协同过滤、HMM 等等