EM（期望最大，Expectation-Maximization）算法是一种常用的迭代优化算法，通常用于含有隐变量或不完全数据的问题中，旨在估计模型的参数，使得对观测数据的对数似然函数达到最大化。它广泛应用于混合高斯模型（GMM）、隐马尔可夫模型（HMM）、协同过滤等问题中。

# EM 算法的基本思想

EM 算法通过迭代地执行两个步骤：

  • E 步（期望步，Expectation Step）：在当前参数的基础上，计算隐含变量的期望值。
  • M 步（最大化步，Maximization Step）：给定隐含变量的期望值，最大化似然函数，重新估计模型参数。

这个过程会在 E 步和 M 步之间反复迭代，直到模型的参数收敛到一个局部最优解。

# EM 算法的核心步骤

假设我们有一些带有隐变量的数据，数据的联合分布为 P (X,Z∣θ)，其中：

  • X 是观测数据（可见数据）。
  • Z 是隐变量（隐藏数据）。
  • θ 是模型的参数，我们希望通过 EM 算法来估计这些参数。

EM 算法的目标是通过最大化对观测数据的似然函数来估计参数：

L(θ)=P(X∣θ)=Z∑​P(X,Z∣θ)L(θ)=P(X∣θ)=Z∑​P(X,Z∣θ)

由于隐变量 Z 的存在，直接求解对数似然比较复杂，因此通过 EM 算法来迭代优化。具体步骤如下：

# 1. 初始化：

首先，我们对参数 θ 进行初始化。这些初始值可以随机选取或通过某些启发式方法得到。

# 2. E 步（期望步）：

在 E 步中，给定当前模型参数 θ(t)，计算隐变量的条件期望。直观上讲，这一步是计算在当前参数下，隐变量 Z 的可能取值。具体而言，E 步计算的是后验分布：

Q(θ∣θ(t))=EZ∣X,θ(t)​[logP(X,Z∣θ)]Q(θ∣θ^{(t)})=E_{Z∣X,θ^{(t)}}​[logP(X,Z∣θ)]

这相当于在当前参数 θ(t) 下，给定观测数据 X，计算 Z 的期望值，从而得到一个新的目标函数 Q (θ) 来表示。

# 3. M 步（最大化步）：

在 M 步中，最大化步骤是通过更新参数 θ 来最大化上一步中计算得到的 Q (θ∣θ(t))：

θ(t+1)=argθmax​Q(θ∣θ(t))θ(t+1)=argθmax​Q(θ∣θ(t))

这一优化过程可以被视为是寻找使得期望似然函数最大的参数。

# 4. 迭代过程：

重复执行 E 步和 M 步，直到参数 θ 的变化足够小（即达到某个收敛条件），或者对数似然函数的改进趋于停止。

# EM 算法的直观理解

  1. 初始参数：在初始参数下，模型对数据的描述可能不太准确，但给定参数后，我们可以计算隐含变量的期望值（E 步）。
  2. 隐变量的更新：基于当前参数的模型假设，推断每个数据点的隐变量（如在高斯混合模型中，这是每个数据点属于哪个高斯成分的概率）。
  3. 参数的更新：在给定隐变量的期望值后，最大化似然估计，更新模型参数，使模型更好地拟合观测数据（M 步）。
  4. 重复：这一过程反复执行，随着参数的更新，模型逐渐收敛到一个能解释数据的较优状态。

# EM 算法的应用示例

高斯混合模型（GMM） 是 EM 算法的典型应用场景。假设数据是由多个高斯分布混合生成的，我们不知道每个数据点来自哪个高斯分布。EM 算法帮助我们估计每个数据点属于哪个高斯分布（隐含变量），以及每个高斯分布的参数。

GMM 中的 EM 算法过程：

  1. 初始化：初始化每个高斯成分的均值、协方差矩阵和混合系数。
  2. E 步：根据当前的高斯成分参数，计算每个数据点属于每个高斯分布的概率（即计算隐含变量的期望值）。
  3. M 步：利用这些概率来更新每个高斯分布的参数（均值、协方差和混合系数）。
  4. 迭代：不断重复 E 步和 M 步，直到参数收敛。

# 优缺点

  • 优点：

    • 适用于隐变量模型，尤其是观测数据不完整或带有噪声的情况。
    • 每次更新步骤都会增加对数似然函数的值，保证算法的收敛性。

  • 缺点：

    • EM 算法只能保证局部最优解，无法保证找到全局最优解。
    • 对初始参数较为敏感，容易陷入局部最优。
    • 在某些模型中，E 步可能计算复杂或涉及数值不稳定性。

# 总结

EM 算法通过两个步骤交替迭代：E 步估计隐变量的期望，M 步最大化期望似然函数，逐步优化模型参数。它是处理隐变量或不完全数据的一种常用方法，广泛应用于聚类分析、混合模型和隐马尔可夫模型等领域。

# 代码实现

# 高斯混合模型（GMM）中的 EM 算法

在这个例子中，我们有两组高斯分布的数据，并希望使用 EM 算法估计它们的参数。

# 代码说明

  1. 数据生成：我们生成了两组二维高斯分布的数据，分别表示两个类。
  2. EM 算法：
    • E 步：在 e_step 函数中，给定当前参数（均值、协方差和权重），计算每个数据点属于每个高斯成分的后验概率（即责任度 responsibilities ）。
    • M 步：在 m_step 函数中，使用 E 步计算的责任度来更新每个高斯成分的参数（包括均值、协方差和混合系数）。
    • 对数似然函数：在 compute_log_likelihood 中，我们计算了当前参数下数据的对数似然，用于判断算法是否收敛。
  3. 迭代过程：EM 算法会不断执行 E 步和 M 步，直到对数似然函数收敛。

# 结果解释

最终输出的参数包括：

  • 混合系数（weights）：表示每个高斯成分的权重，即每个类的比例。
  • 均值（means）：每个高斯分布的中心位置。
  • 协方差矩阵（covariances）：描述每个高斯分布的形状和大小。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
import numpy as np
from scipy.stats import multivariate_normal

# 随机生成两个高斯分布的数据
np.random.seed(0)
n_samples = 300

# 生成数据
mean1 = [0, 0]
cov1 = [[1, 0], [0, 1]]  # 协方差矩阵为单位矩阵
data1 = np.random.multivariate_normal(mean1, cov1, int(0.5 * n_samples))

mean2 = [3, 3]
cov2 = [[1, 0.2], [0.2, 1]]  # 协方差矩阵有一定的相关性
data2 = np.random.multivariate_normal(mean2, cov2, int(0.5 * n_samples))

# 合并两组数据
data = np.vstack((data1, data2))

# EM 算法参数初始化
def initialize_parameters(n_clusters, n_features):
    weights = np.ones(n_clusters) / n_clusters  # 初始化每个高斯分布的权重
    means = np.random.rand(n_clusters, n_features) * 5  # 初始化均值
    covariances = np.array([np.eye(n_features)] * n_clusters)  # 初始化协方差矩阵为单位矩阵
    return weights, means, covariances

# 计算 E 步，返回每个点属于每个高斯成分的后验概率
def e_step(data, weights, means, covariances, n_clusters):
    n_samples = data.shape[0]
    responsibilities = np.zeros((n_samples, n_clusters))

    for k in range(n_clusters):
        rv = multivariate_normal(mean=means[k], cov=covariances[k])
        responsibilities[:, k] = weights[k] * rv.pdf(data)
    
    # 归一化 responsibilities，使每个点属于不同高斯分布的概率和为1
    responsibilities /= responsibilities.sum(axis=1, keepdims=True)
    return responsibilities

# 计算 M 步，更新模型参数（均值、协方差和混合系数）
def m_step(data, responsibilities, n_clusters):
    n_samples, n_features = data.shape
    weights = np.sum(responsibilities, axis=0) / n_samples
    means = np.dot(responsibilities.T, data) / np.sum(responsibilities, axis=0)[:, np.newaxis]
    
    covariances = np.zeros((n_clusters, n_features, n_features))
    for k in range(n_clusters):
        diff = data - means[k]
        covariances[k] = np.dot(responsibilities[:, k] * diff.T, diff) / np.sum(responsibilities[:, k])
    
    return weights, means, covariances

# 计算对数似然函数，判断收敛
def compute_log_likelihood(data, weights, means, covariances, n_clusters):
    log_likelihood = 0
    for k in range(n_clusters):
        rv = multivariate_normal(mean=means[k], cov=covariances[k])
        log_likelihood += weights[k] * rv.pdf(data)
    return np.log(log_likelihood).sum()

# EM 算法主函数
def em_algorithm(data, n_clusters, n_iter=100, tol=1e-4):
    n_samples, n_features = data.shape
    weights, means, covariances = initialize_parameters(n_clusters, n_features)
    
    log_likelihoods = []
    
    for i in range(n_iter):
        # E 步：计算 responsibilities
        responsibilities = e_step(data, weights, means, covariances, n_clusters)
        
        # M 步：更新权重、均值和协方差
        weights, means, covariances = m_step(data, responsibilities, n_clusters)
        
        # 计算对数似然
        log_likelihood = compute_log_likelihood(data, weights, means, covariances, n_clusters)
        log_likelihoods.append(log_likelihood)
        
        # 检查收敛
        if i > 0 and np.abs(log_likelihoods[-1] - log_likelihoods[-2]) < tol:
            print(f"EM算法在第 {i} 次迭代时收敛。")
            break
    
    return weights, means, covariances, responsibilities

# 运行EM算法，寻找高斯混合模型的参数
n_clusters = 2
weights, means, covariances, responsibilities = em_algorithm(data, n_clusters)

print("混合系数:", weights)
print("均值:", means)
print("协方差:", covariances)