EM算法到底是什么东东

frostmelody

已于 2025-04-06 11:41:24 修改

阅读量1.2k

点赞数 18

CC 4.0 BY-SA版权

分类专栏：机器学习小知识点 PyTorch小知识点文章标签：算法机器学习人工智能

于 2025-04-06 01:03:32 首次发布

本文链接：https://blog.youkuaiyun.com/Listennnn/article/details/147016470

机器学习小知识点同时被 2 个专栏收录

16 篇文章

订阅专栏

PyTorch小知识点

16 篇文章

订阅专栏

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

🐤 第一步：直观理解

EM算法的核心是：

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

例子1：

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

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

EM算法的做法是：

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

✍️ 第二步：数学公式

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

每个分布 $k$ 有自己的参数：均值 $μk\mu_k$ 、方差 $σk2\sigma_k^2$ 、权重 $πk\pi_k$ （概率总和为1）
但你不知道哪个点来自哪个分布（这是隐变量）

E步（Expectation），即“先猜”

初始化：随机初始化均值 $μk\mu_k$ 、方差 $σk2\sigma_k^2$ 和权重 $πk\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），即“反推参数”

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

$\mu_k = \frac{\sum_i \gamma_{ik} x_i}{\sum_i \gamma_{ik}}, \quad \sigma_k^2 = \frac{\sum_i \gamma_{ik} (x_i - \mu_k)^2}{\sum_i \gamma_{ik}}, \quad \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等等，但我们不知道每个学生是小学生还是中学生。我们猜测这些身高来自两个群体：

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

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

弄清楚每个学生属于小学生还是中学生的概率。
估计两个群体的参数：比例（ $π$ ）、均值（ $μ$ ）、标准差（ $σ$ ）。

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

EM算法是什么？

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

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

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

具体案例：一步步拆解

1. 初始化：随便猜参数

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

小学生：
- 比例（ $π_1$ ）：50%（0.5）
- 均值（ $μ_1$ ）：150cm
- 标准差（ $σ_1$ ）：5cm
中学生：
- 比例（ $π_2$ ）：50%（0.5）
- 均值（ $μ_2$ ）：170cm
- 标准差（ $σ_2$ ）：6cm

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

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

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

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

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

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

正态分布的公式是：
$P(X)=12π⋅σexp⁡(−(X−μ)22σ2)P(X)=\frac{1}{\sqrt{2\pi}\cdot\sigma}\exp\left(-\frac{(X-\mu)^2}{2\sigma^2}\right)$

小学生：
$P(小学生∣160)=12π⋅5exp⁡(−(160−150)22⋅52)P(\text{小学生}|160)=\frac{1}{\sqrt{2\pi}\cdot5}\exp\left(-\frac{(160-150)^2}{2\cdot5^2}\right)$
计算指数部分：(160 - 150)² = 100，2 × 5² = 50，-100 / 50 = -2，exp(-2) ≈ 0.135。所以结果是一个较小的数。
中学生：
$P(中学生∣160)=12π⋅6exp⁡(−(160−170)22⋅62)P(\text{中学生}|160)=\frac{1}{\sqrt{2\pi}\cdot6}\exp\left(-\frac{(160-170)^2}{2\cdot6^2}\right)$
计算指数部分：(160 - 170)² = 100，2 × 6² = 72，-100 / 72 ≈ -1.39，exp(-1.39) ≈ 0.25。结果比小学生的稍大。

简单来说：

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

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

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

假设计算后：

$P(小学生∣160)≈0.3P(\text{小学生}|160)≈0.3$ （30%）
$P(中学生∣160)≈0.7P(\text{中学生}|160)≈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）更新比例（ $π$ ）

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

（2）更新均值（ $μ$ ）

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

（3）更新标准差（ $σ$ ）

新 $σ_1=$ 身高偏离新均值 $μ_1$ 的加权方差：
$σ1=0.3⋅(160−158.5)2+0.6⋅(155−158.5)2+0.1⋅(175−158.5)21.0σ_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。
新 $σ_2=$ 类似计算，得到新值，比如5.5cm。

4. 重复迭代

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

小学生： $π1≈0.4,μ1≈148cm,σ1≈4cmπ_1≈0.4,μ_1≈148\text{cm},σ_1≈4\text{cm}$
中学生： $π2≈0.6,μ2≈172cm,σ2≈5cmπ_2≈0.6,μ_2≈172\text{cm},σ_2≈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 plt


class 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=100
            EM算法最大迭代次数
        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.0
            for 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
        # 归一化使各学生概率和为1
        responsibilities /= responsibilities.sum(axis=1, keepdims=True)
        return responsibilities

    def _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) ** 2
        self.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 in
                   zip([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 = False

        for 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 = True
                break

        # 最终按身高均值排序群体（保证小学生群体在前）
        order = np.argsort(self.mu)
        self.mu = self.mu[order]
        self.pi = self.pi[order]
        self.sigma = self.sigma[order]

        return self

    def 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 = model
        self.data = data
        self.true_labels = true_labels
        self.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 == k
            print(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 等等