EM算法与高斯混合模型（GMM）系统梳理

一、核心前置概念

在学习EM和GMM前，需先明确以下基础概念，是理解后续内容的关键：

概念	定义与作用
隐变量（Latent Variable）	无法直接观测但影响可观测变量的变量，是EM算法的核心解决对象（如GMM中样本所属的高斯成分）
最大似然估计（MLE）	给定观测数据，求解使似然函数最大的模型参数，适用于无隐变量的场景
高斯分布（正态分布）	单变量：$\mathcal{N}(x\mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$；多变量：$\mathcal{N}(\mathbf{x}\mathbf{\mu },\mathbf{\Sigma })=\frac{1}{(2\pi {)}^{d/2}{\mathbf{\Sigma }}^{1/2}}{e}^{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })}$（$d$为特征维度）

二、EM算法（期望最大化算法）

1. 算法定义

EM（Expectation-Maximization）是针对含隐变量的概率模型的参数估计迭代算法，核心是通过交替执行“期望步（E步）”和“最大化步（M步）”，逼近含隐变量模型的最大似然解（MLE）。

2. 核心思想

含隐变量的对数似然函数$\log p(X|\Theta )$（$X$为观测数据，$\Theta$为模型参数）难以直接最大化，EM通过构造证据下界（ELBO） 替代原目标，交替优化下界以逼近最优解，保证每轮迭代后似然值非递减。

3. 通用迭代步骤

步骤	操作
初始化	设定模型参数初始值${\Theta }^{(0)}$（如随机初始化），设定收敛阈值$ϵ$
E步	计算隐变量的后验概率，构造对数似然关于隐变量的期望（Q函数）： $Q(\Theta {\Theta }^{(t)})={\mathbb{E}}_{ZX,{\Theta }^{(t)}}[\log p(X,Z\Theta )]$
M步	最大化Q函数，更新参数：${\Theta }^{(t+1)}=\arg {\max }_{\Theta }Q(\Theta {\Theta }^{(t)}$
收敛判断	若$|{\Theta }^{(t+1)}-{\Theta }^{(t)}|<ϵ$或似然值变化<$ϵ$，停止迭代；否则回到E步

4. 收敛性

EM算法保证每轮迭代后对数似然值非递减，最终收敛到局部最优解（而非全局最优，初始值会影响结果）。

三、高斯混合模型（GMM）

1. 模型定义

GMM是由$K$个高斯分布（称为“成分”）加权组合而成的概率模型，用于建模多峰、非高斯的复杂数据分布，公式如下：

（1）单变量GMM

$$p(x)=\sum _{k=1}^K{\pi }_k\mathcal{N}(x|{\mu }_k,{\sigma }_k^2)$$

（2）多变量GMM

$$p(\mathbf{x})=\sum _{k=1}^K{\pi }_k\mathcal{N}(\mathbf{x}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)$$

其中核心约束与参数说明：

• ${\pi }_k$：混合系数（第$k$个成分的权重），满足${\sum }_{k=1}^K{\pi }_k=1$且${\pi }_k\ge 0$；
• ${\mathbf{\mu }}_k$：第$k$个成分的$d$维均值向量；
• ${\mathbf{\Sigma }}_k$：第$k$个成分的$d\times d$协方差矩阵（正定）；
• 模型总参数：Θ={π1,...,πK,μ1,...,μK,Σ1,...,ΣK}\Theta = \{\pi_1,...,\pi_K, \boldsymbol{\mu}_1,...,\boldsymbol{\mu}_K, \boldsymbol{\Sigma}_1,...,\boldsymbol{\Sigma}_K\}Θ={π1​,...,πK​,μ1​,...,μK​,Σ1​,...,ΣK​}。

2. 用EM算法求解GMM（核心步骤）

GMM的参数估计是EM算法的经典应用，核心是将“样本所属的高斯成分”作为隐变量$z_{nk}$（$z_{nk}=1$表示样本$n$属于第$k$个成分，${\sum }_k{z}_{nk}=1$）。

（1）E步：计算“责任（Responsibility）”

责任${\gamma }_{nk}$表示给定当前参数和样本$x_n$时，样本$n$属于第$k$个成分的后验概率：
$${\gamma }_{nk}=p(z_{nk}=1|x_n,{\Theta }^{(t)})=\frac{{\pi }_k^{(t)}\mathcal{N}(x_n|{\mathbf{\mu }}_k^{(t)},{\mathbf{\Sigma }}_k^{(t)})}{{\sum }_{j=1}^K{\pi }_j^{(t)}\mathcal{N}(x_n|{\mathbf{\mu }}_j^{(t)},{\mathbf{\Sigma }}_j^{(t)})}$$

（2）M步：更新模型参数

基于E步的责任，最大化Q函数，得到参数更新公式：

参数类型	更新公式
混合系数	${\pi }_k^{(t+1)}=\frac{1}{N}{\sum }_{n=1}^N{\gamma }_{nk}$（$N$为样本总数）
均值	${\mathbf{\mu }}_k^{(t+1)}=\frac{{\sum }_{n=1}^N{\gamma }_{nk}{\mathbf{x}}_n}{{\sum }_{n=1}^N{\gamma }_{nk}}$
协方差	${\mathbf{\Sigma }}_k^{(t+1)}=\frac{{\sum }_{n=1}^N{\gamma }_{nk}({\mathbf{x}}_n-{\mathbf{\mu }}_k^{(t+1)})({\mathbf{x}}_n-{\mathbf{\mu }}_k^{(t+1)}{)}^T}{{\sum }_{n=1}^N{\gamma }_{nk}}$

（3）收敛条件

当参数变化的L2范数小于阈值（如$10^{-6}$），或对数似然值$\log p(X|\Theta )$的变化小于阈值时，停止迭代。

3. GMM vs K-Means（关键对比）

GMM常与K-Means（硬聚类）对比，核心差异如下：

维度	GMM（高斯混合模型）	K-Means（K均值聚类）
聚类类型	软聚类（输出样本属于每个成分的概率）	硬聚类（输出样本唯一所属类别）
模型假设	数据服从高斯混合分布，考虑均值+协方差	数据簇为球形，仅考虑均值（距离）
求解方法	EM算法	迭代重分配（距离最小化）
输出结果	概率分布+聚类标签	仅聚类标签
鲁棒性	对异常值更敏感（协方差易受影响）	对异常值敏感（距离易受影响）

四、Python实现语法与案例

1. 核心API语法（sklearn）

sklearn提供GaussianMixture类实现GMM，核心参数与方法如下：

类别	内容
导入	from sklearn.mixture import GaussianMixture
初始化	gmm = GaussianMixture(n_components=K, # 高斯成分数<br> covariance_type='full', # 协方差类型：full/ diag/ spherical/ tied<br> max_iter=100, # 最大迭代次数<br> tol=1e-6, # 收敛阈值<br> random_state=42)
训练	gmm.fit(X) # X为形状(n_samples, n_features)的数组
预测聚类标签	labels = gmm.predict(X) # 硬聚类标签（概率最大的成分）
预测概率	probs = gmm.predict_proba(X) # 形状(n_samples, K)，每个样本属于各成分的概率
模型参数	gmm.weights_ # 混合系数π gmm.means_ # 均值μ gmm.covariances_ # 协方差Σ

2. 完整案例：GMM拟合二维模拟数据

步骤1：生成模拟数据

import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
from scipy.stats import multivariate_normal

# 1. 生成二维高斯混合数据
np.random.seed(42)
# 定义3个高斯成分的参数
K = 3
n_samples = 500
# 成分1：均值[0,0]，协方差[[1,0.5],[0.5,1]]，样本数200
mean1 = [0, 0]
cov1 = [[1, 0.5], [0.5, 1]]
data1 = np.random.multivariate_normal(mean1, cov1, 200)
# 成分2：均值[5,5]，协方差[[1,-0.5],[-0.5,1]]，样本数150
mean2 = [5, 5]
cov2 = [[1, -0.5], [-0.5, 1]]
data2 = np.random.multivariate_normal(mean2, cov2, 150)
# 成分3：均值[0,8]，协方差[[1,0],[0,1]]，样本数150
mean3 = [0, 8]
cov3 = [[1, 0], [0, 1]]
data3 = np.random.multivariate_normal(mean3, cov3, 150)
# 合并数据
X = np.vstack([data1, data2, data3])
y_true = np.hstack([np.zeros(200), np.ones(150), np.ones(150)*2])  # 真实标签

# 可视化原始数据
plt.figure(figsize=(8, 6))
plt.scatter(X[:,0], X[:,1], c=y_true, cmap='viridis', alpha=0.6)
plt.title('原始模拟数据（3个高斯成分）')
plt.xlabel('特征1')
plt.ylabel('特征2')
plt.show()

步骤2：sklearn实现GMM

# 2. 用sklearn拟合GMM
gmm = GaussianMixture(
    n_components=3,  # 3个成分
    covariance_type='full',  # 全协方差矩阵
    max_iter=200,
    tol=1e-6,
    random_state=42
)
gmm.fit(X)

# 输出模型参数
print("=== GMM模型参数 ===")
print(f"混合系数（π）：{gmm.weights_}")
print(f"均值（μ）：\n{gmm.means_}")
print(f"协方差（Σ）：\n{gmm.covariances_}")

# 预测聚类标签和概率
labels_pred = gmm.predict(X)
probs_pred = gmm.predict_proba(X)

# 可视化拟合结果
plt.figure(figsize=(8, 6))
# 绘制样本
plt.scatter(X[:,0], X[:,1], c=labels_pred, cmap='viridis', alpha=0.6)
# 绘制每个高斯成分的等高线
x_min, x_max = X[:,0].min()-1, X[:,0].max()+1
y_min, y_max = X[:,1].min()-1, X[:,1].max()+1
xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100), np.linspace(y_min, y_max, 100))
zz = np.zeros_like(xx)
for k in range(K):
    # 计算每个网格点的概率密度
    rv = multivariate_normal(gmm.means_[k], gmm.covariances_[k])
    zz += gmm.weights_[k] * rv.pdf(np.dstack((xx, yy)))
# 绘制等高线
plt.contour(xx, yy, zz, levels=5, colors='black', alpha=0.5)
plt.title('GMM拟合结果（sklearn）')
plt.xlabel('特征1')
plt.ylabel('特征2')
plt.show()

步骤3：手动实现EM算法拟合GMM（核心逻辑）

def gmm_em(X, K, max_iter=200, tol=1e-6):
    """
    手动实现EM算法拟合GMM
    参数：
        X: 数据，形状(n_samples, n_features)
        K: 高斯成分数
        max_iter: 最大迭代次数
        tol: 收敛阈值
    返回：
        pi: 混合系数，形状(K,)
        mu: 均值，形状(K, n_features)
        sigma: 协方差，形状(K, n_features, n_features)
        gamma: 责任，形状(n_samples, K)
    """
    n_samples, n_features = X.shape
    
    # 1. 初始化参数
    # 混合系数：均匀初始化
    pi = np.ones(K) / K
    # 均值：随机选择K个样本作为初始均值
    idx = np.random.choice(n_samples, K, replace=False)
    mu = X[idx]
    # 协方差：初始化为单位矩阵
    sigma = np.array([np.eye(n_features) for _ in range(K)])
    
    # 初始化对数似然
    log_likelihood = []
    
    for _ in range(max_iter):
        # 2. E步：计算责任gamma
        gamma = np.zeros((n_samples, K))
        for k in range(K):
            # 计算每个样本在第k个成分的概率密度
            gamma[:, k] = pi[k] * multivariate_normal.pdf(X, mean=mu[k], cov=sigma[k])
        # 归一化得到后验概率
        gamma = gamma / gamma.sum(axis=1, keepdims=True)
        
        # 计算当前对数似然
        ll = np.sum(np.log(np.sum([pi[k] * multivariate_normal.pdf(X, mu[k], sigma[k]) for k in range(K)], axis=0)))
        log_likelihood.append(ll)
        
        # 3. M步：更新参数
        # 更新混合系数
        pi_new = gamma.sum(axis=0) / n_samples
        # 更新均值
        mu_new = np.zeros_like(mu)
        for k in range(K):
            mu_new[k] = (gamma[:, k].reshape(-1,1) * X).sum(axis=0) / gamma[:, k].sum()
        # 更新协方差
        sigma_new = np.zeros_like(sigma)
        for k in range(K):
            diff = X - mu_new[k]
            sigma_new[k] = (gamma[:, k].reshape(-1,1,1) * np.einsum('ni,nj->nij', diff, diff)).sum(axis=0) / gamma[:, k].sum()
        
        # 4. 收敛判断
        if np.linalg.norm(pi_new - pi) < tol and \
           np.linalg.norm(mu_new - mu) < tol and \
           np.linalg.norm(sigma_new - sigma) < tol:
            break
        
        # 更新参数
        pi, mu, sigma = pi_new, mu_new, sigma_new
    
    return pi, mu, sigma, gamma

# 调用手动实现的EM-GMM
pi_manual, mu_manual, sigma_manual, gamma_manual = gmm_em(X, K=3, max_iter=200, tol=1e-6)

# 输出手动拟合的参数
print("\n=== 手动EM-GMM模型参数 ===")
print(f"混合系数（π）：{pi_manual}")
print(f"均值（μ）：\n{mu_manual}")
print(f"协方差（Σ）：\n{sigma_manual}")

# 可视化手动拟合结果
labels_manual = np.argmax(gamma_manual, axis=1)
plt.figure(figsize=(8, 6))
plt.scatter(X[:,0], X[:,1], c=labels_manual, cmap='viridis', alpha=0.6)
# 绘制等高线
x_min, x_max = X[:,0].min()-1, X[:,0].max()+1
y_min, y_max = X[:,1].min()-1, X[:,1].max()+1
xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100), np.linspace(y_min, y_max, 100))
zz = np.zeros_like(xx)
for k in range(K):
    rv = multivariate_normal(mu_manual[k], sigma_manual[k])
    zz += pi_manual[k] * rv.pdf(np.dstack((xx, yy)))
plt.contour(xx, yy, zz, levels=5, colors='black', alpha=0.5)
plt.title('手动EM算法拟合GMM结果')
plt.xlabel('特征1')
plt.ylabel('特征2')
plt.show()

3. 案例结果说明

• 模拟数据由3个二维高斯成分组成，sklearn和手动实现的EM-GMM均能较好拟合出3个成分的参数；
• 手动实现的参数与sklearn结果略有差异（初始值不同），但聚类效果一致；
• 等高线清晰展示了每个高斯成分的分布范围，混合系数反映了各成分的样本占比。

五、关键注意事项

1. 成分数K的选择：需通过信息准则（如AIC、BIC）确定，sklearn的GaussianMixture支持aic()和bic()方法；
2. 初始值敏感：EM易陷入局部最优，可多次随机初始化取似然值最大的结果；
3. 协方差类型：covariance_type参数（full/diag/spherical/tied）需根据数据特性选择（如球形协方差等价于K-Means的球形假设）；
4. 数据预处理：GMM对数据尺度敏感，建议先标准化（StandardScaler）。

六、总结

EM算法是处理含隐变量模型参数估计的通用框架，GMM是EM算法的经典应用，核心是通过“责任”的交替计算与参数更新拟合多峰分布。相比K-Means，GMM能输出概率分布，更贴合真实数据的统计特性，是机器学习中聚类、密度估计、生成模型的重要工具。