EM算法——含有隐变量的最大似然估计

EM算法（Expectation-Maximization）是一种用于解决含有隐变量的概率模型参数估计问题的迭代优化算法。其基本思想是通过交替进行期望（Expectation）和最大化（Maximization）两个步骤来优化模型参数。在E步骤中，通过当前参数对隐变量的条件分布进行估计，计算完全数据对数似然的期望值。这一步旨在估计隐变量的信息。在M步骤中，利用E步骤中得到的隐变量的期望值最大化对数似然函数，从而更新模型参数。这两个步骤交替进行直至收敛。EM算法广泛应用于许多领域，其中最典型的包括聚类、混合模型、隐马尔可夫模型等。在聚类中，EM算法可用于估计混合高斯模型的参数，将每个观测数据点分配到各个高斯分布中。在隐马尔可夫模型中，EM算法可以用于学习模型的转移概率和发射概率。EM算法的优势在于能够处理含有不完全数据或缺失数据的情况，使其在实际问题中得到广泛应用。其迭代的过程不断提高对未观测变量的估计精度，从而提高对模型参数的估计准确性，为复杂问题提供了一种有效的解决方案。

一、EM算法引入

假定有两枚不同的硬币A和B，它们的重量分布（跑一次出现正面概率） $\theta_A$ 和 $\theta_B$ 是未知的，其数值可通过抛掷后计算正反面各自出现的次数来估计。具体的估计方法是在每一轮中随机抽出一枚硬币抛掷10次，同样的过程执行5轮，根据这50次投币的结果来计算 $\theta_A$ 和 $\theta_B$ 的最大似然估计。

图1	图2

在上图的单次试验中，硬币A被抽到3次，30次投掷中出现了24次正面；硬币B被抽到2次，20次投掷中出现了9次正面。用最大似然估计可以计算出

{\hat{θ}}_{A} = 24 / (24 + 6) = 0.8, {\hat{θ}}_{B} = 9 / (9 + 11) = 0.45

$\hat{\theta}_A=24/(24+6)=0.8,\hat{\theta}_B=9/(9+11)=0.45$

这个问题很简单，但是如果我们只能知道每一轮中出现的正反面结果，却不能得知到底选取的硬币是A还是B，问题就没那么简单了。这里的硬币选择就是不能直接观测的隐变量。如果不管这个隐变量，就无法估计未知参数；要确定这一组隐变量，又得基于未知的硬币重量分布，问题进入了死胡同。既然数据中信息不完整，那就人为地补充完整。在这个问题中，隐藏的硬币选择和待估计的重量分布，两者确定一个就可以计算另一个。
由于观测结果给出了关于重量分布的信息，那就不放人为设定一组初始化的参数 $\hat{\theta}^{(t)}=( \hat{\theta}_A^{(t)},\hat{\theta}_B^{(t)} )$ ，用这组猜测的重量分布去倒推到底每一轮使用的是哪个硬币。计算出来的硬币选择会被用来对原来随机产生的初始化参数进行更新。如果硬币选择的结果是正确的，就可以利用最大似然估计计算出新的参数 $\hat{\theta}^{(t+1)}$ 。而更新后的参数有可以应用在观测结果上，对硬币选择的结果进行修正，从而形成了“批评-自我批评”的循环过程。这个过程会持续到隐变量和未知参数的取值都不再发生变化，其结果就是最终的输出。
将思路应用到上图的投掷结果中，就是EM算法的雏形。若两个初始的参数随机设定为 $\hat{\theta}_A^{(0)}=0.6,\hat{\theta}_B^{(0)}=0.5$ ，在这两个参数下出现第一轮结果，也就是5正5反的概率就可以表示为

P (H^{5} T^{5} | A) = {0.6}^{5} \times {0.4}^{5}, P (H^{5} T^{5} | B) = {0.5}^{10}

$P(H^5T^5|A)=0.6^5 \times 0.4^5,P(H^5T^5|B)=0.5^{10}$

对上面的两个似然概率进行归一化可以得出后验概率，两者分别是0.45和0.55，也就是上图中的结果。这说明如果初始结果的随机参数是准确的，那第一轮结果更可能由硬币B产生（0.55>0.45）。同理也可以计算出其他四轮的结果来自不同硬币的后验概率，结果已经在上图中显示。
在已知硬币的选择时，所有正反面的结果都有明确的归属：要么来自A要么来自B。利用后验概率可以直接对硬币的选择作出判断：1、4轮使用的是硬币B，2、3、5轮使用的是硬币A。既然硬币的选择已经确定，这时就可以使用最大似然估计，其结果和前文中的最大似然估计结果是相同的，也就是 $\hat{\theta}_A^{(1)}=0.8, \quad \hat{\theta}_B^{(1)}=0.45$ 。利用这组更新的参数又可以重新计算每一轮次抽取不同硬币的后验概率。虽然这种方法能够实现隐变量和参数的动态更新，但它还不是真正的EM算法，而是硬输出的k均值聚类。真正的EM算法并不会将后验概率最大的值赋给隐变量，而是考虑其所有可能的取值，在概率分布的框架下进行分析。
在前面的例子中，由于第一轮投掷硬币A的可能性是0.45，那么硬币A对正反面出现次数的贡献就是45%，在5次正面的结果中，来源于硬币A的就是 5×0.45=2.25 次，来源于硬币B的则是2.75次。同理可以计算出其他轮次中A和B各自的贡献，贡献的比例都和计算出的后验概率相对应。计算出A和B在不同轮次中的贡献，就可以对未知参数做出更加精确的估计。在 50 次投掷中，硬币A贡献了21.3次正面和8.6次反面，其参数估计值 $\hat{\theta}_A^{(1)}=0.71$ ；硬币B贡献了11.7次正面和8.4次反面，其参数估计值 $\hat{\theta}_B^{(1)}=0.58$ 。利用这组参数继续迭代更新，就可以计算出最终的估计值。

二、EM算法

对于上面硬币重量估计，多了一个硬币种类类型的隐变量，设为 $z$ ，可以把它看做一个10维的向量 $(z_1, z_2,...,z_{10})$ ，代表每次投掷时所使用的硬币，比如 $z_1$ ，就代表第一轮投掷时使用的硬币是A还是B。但是，这个变量 $Z$ 不知道，就无法去估计 $P_A$ 和 $P_B$ ，所以我们先必须估计出 $Z$ ，然后才能进一步估计 $P_A$ 和 $P_B$ 。可是要估计 $Z$ ，我们又得知道 $P_A$ 和 $P_B$ ，这样我们才能用极大似然概率法去估计 $Z$ ，这不是前面提的鸡生蛋和蛋生鸡的问题吗，如何解？答案就是先随机初始化一个 $P_A$ 和 $P_B$ ，用它来估计 $Z$ ，然后基于 $Z$ ，还是按照最大似然概率法则去估计新的 $P_A$ 和 $P_B$ ，如果新的 $P_A$ 和 $P_B$ 和我们初始化的 $P_A$ 和 $P_B$ 一样，请问这说明了什么？这说明我们初始化的 $P_A$ 和 $P_B$ 是一个相当靠谱的估计！这就是说，我们初始化的 $P_A$ 和 $P_B$ ，按照最大似然概率就可以估计出 $Z$ ，然后基于 $Z$ ，按照最大似然概率可以反过来估计出 $P_1$ 和 $P_2$ ，当与我们初始化的 $P_A$ 和 $P_B$ 一样时，说明是 $P_1$ 和 $P_2$ 很可能就是真实的值。这里包含了两个交互的最大似然估计。如果新估计出来的 $P_A$ 和 $P_B$ 与我们初始化的值差别很大，怎么办呢？就是继续用新的 $P_1$ 和 $P_2$ 迭代，直至收敛。
概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，我们讨论极大似然估计，极大后验概率估计与其类似。
输入：观测变量数据Y，隐变量数据Z，联合分布 $P(Y,Z|\theta)$ ，条件分布 $P(Z|Y,\theta)$ ；输出：模型参数 $\theta$

(1)选择参数的初值 $\theta^0$ ，开始迭代；
(2)E步：记 $\theta^i$ 为第 $i$ 次迭代参数 $\theta$ 的估计值，在第 $i+1$ 次迭代的E步，计算

\begin{aligned} Q (θ, θ^{i}) & = E_{Z} [l o g P (Y, Z | θ) | Y, θ^{i}] \\ = \sum_{Z} l o g P (Y, Z | θ) P (Z | Y, θ^{i}) \end{aligned}

$\begin{aligned} Q(\theta,\theta^i)&=E_{Z}[logP(Y,Z|\theta)|Y,\theta^i]\\ & =\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^i) \end{aligned}$

这里， $P(Z|Y,\theta^i)$ 是在给定观测数据Y和当前的参数估计 $\theta^i$ 下隐变量数据Z的条件概率分布；

(3) M步：求使 $Q(\theta,\theta^i)$ 极大化的 $\theta$ ，确定第 $i+1$ 次迭代的参数的估计值 $\theta^{i+1}$ ，
$\theta^{i+1}=arg \max \limits_{\theta}Q(\theta,\theta^{i})$

$Q(\theta,\theta^{i})$ 是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的。

(4) 重复第(2)步和第(3)步，直到收敛，收敛条件：
$|| \theta^{i+1}-\theta^{i} || < \varepsilon_1$

或者：

| | Q (θ^{i + 1}, θ^{i}) - Q (θ^{i}, θ^{i}) | | < ε_{2}

$||Q(\theta^{i+1},\theta^{i})-Q(\theta^{i},\theta^{i})|| <\varepsilon_2$

收敛迭代就结束了。

三、在高斯混合分布中的应用

EM算法的一个重要应用场景就是高斯混合模型的参数估计。

3.1高斯混合分布

现实采集的数据是比较复杂的，通常无法只用一个高斯分布拟合，而是可以看作多个随机过程的混合。可定义高斯混合模型是 $K$ 个高斯分布的组合，用以拟合复杂数据。假设有一个数据集，包含了 $N$ 个相互独立的数据 $x = {x_1, x_2 ...x_i... x_N}$ , 这些数据看起来有 $K$ 个峰，这样的数据集可用以下定义的高斯混合模型拟合：

p (x | Θ) = \sum_{k}^{} α_{k} N (x; μ_{k}, σ_{k}) = \sum_{k}^{} α_{k} \frac{1}{\sqrt{2 π σ_{k}^{2}}} \exp [- \frac{(x - μ_{k})^{2}}{2 σ_{k}^{2}}]

$p(x|\Theta) = \sum_{k}^{}\alpha_{k}\textit{N}(x; \mu_k, \sigma_k) = \sum_{k}^{}\alpha_{k}\frac{1}{\sqrt{2\pi\sigma_k^2}}\exp\left[-\frac{(x-\mu_k)^2}{2\sigma_k^2}\right]$

如果每一个数据点 $x_i$ 都是d维的, 这些数据 $x$ 如上图看起来分散在 $K$ 个聚类，这种数据集可以用多变量高斯混合模型拟合。

p (x | Θ) = \sum_{k}^{} α_{k} N (x; μ_{k}, Σ_{k}) \sum_{k}^{} α_{k} \frac{1}{\sqrt{(2 π)^{d} det (Σ_{k})}} \exp [- \frac{1}{2} (x - μ_{k})^{T} Σ_{k}^{- 1} (x - μ_{k})]

$p(x|\Theta) = \sum_{k}^{}\alpha_{k}\textit{N}(x; \mu_k, \Sigma_k) \nonumber \sum_{k}^{}\alpha_{k}\frac{1}{\sqrt{(2\pi)^d\det(\Sigma_k)}}\exp\left[-\frac{1}{2}(x-\mu_k)^T\Sigma_{k}^{-1}(x-\mu_k)\right]$

其中 $\Theta$ 代表全体高斯模型参数, $\alpha_k$ 是第 $k$ 个高斯模型的先验概率, 各个高斯模型的先验概率加起来等于1，亦即各个高斯分布所占的权重。

\sum_{k} α_{k} = 1

$\sum_{k}\alpha_k = 1$

3.2 EM在高斯混合分布的应用

高斯混合模型就是由多个高斯模型组合在一起的混合模型（可以理解为多个高斯分布函数的线性组合，理论上高斯混合模型是可以拟合任意类型的分布），例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的：

图3	图4

考虑多个高斯模型的混合分布，其中每个高斯模型的权重为 $a_k$ ，均值为 $\mu_k$ ，方差为 $\delta_k^2$ 。观测数据 $y_j$ 可能来自不同的高斯模型，对应的隐藏数据为 $\gamma_{jk}$ ，表示观测数据 $y_j$ 来自第 $k$ 个模型的概率。

初始化:初始化每个高斯模型的参数：权重 $a_k$ 、均值 $\mu_k$ 、方差 $\delta_k^2$ 。
E步骤（Expectation Step）:计算每个观测数据 $y_j$ 来自每个高斯模型的后验概率 $\hat{\gamma_{jk}}$ ：

\hat{γ_{j k}} = \frac{a_{k} ϕ (y_{j} | μ_{k}, δ_{k})}{\sum_{k = 1}^{K} a_{k} ϕ (y_{j} | μ_{k}, δ_{k})}

$\hat{\gamma_{jk}} = \frac{a_k \phi(y_j | \mu_k, \delta_k)}{\sum_{k=1}^{K} a_k \phi(y_j | \mu_k, \delta_k)}$

其中 $\phi(y_j | \mu_k, \delta_k)$ 表示第 $k$ 个高斯分布的密度函数。

M步骤（Maximization Step）:使用E步骤中计算的后验概率 $\hat{\gamma_{jk}}$ 更新高斯模型的参数：
使用E步骤中计算的后验概率 $\hat{\gamma_{jk}}$ 更新高斯模型的参数：

\hat{m_{k}} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}} y_{j}}{\sum_{j = 1}^{N} \hat{γ_{j k}}} {\hat{δ_{k}}}^{2} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}} (y_{j} - \hat{m_{k}})^{2}}{\sum_{j = 1}^{N} \hat{γ_{j k}}} \hat{a_{k}} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}}}{N}

$\hat{m_k} = \frac{\sum_{j=1}^{N} \hat{\gamma_{jk}} y_j}{\sum_{j=1}^{N} \hat{\gamma_{jk}}}\\ \hat{\delta_k}^2 = \frac{\sum_{j=1}^{N} \hat{\gamma_{jk}}(y_j - \hat{m_k})^2}{\sum_{j=1}^{N} \hat{\gamma_{jk}}} \\\hat{a_k} = \frac{\sum_{j=1}^{N} \hat{\gamma_{jk}}}{N}$

更新模型参数:更新模型参数

θ ： θ^{i + 1} = \arg max_{θ} Q (θ, θ^{i})

$\theta：\theta^{i+1} = \arg \max_{\theta} Q(\theta, \theta^i)$

其中 $Q(\theta, \theta^i)$ 是完全似然函数的期望。

收敛判断:判断算法是否收敛，可以通过观察参数的变化或设置收敛阈值。
重复E步骤和M步骤:重复执行E步骤和M步骤直到满足收敛条件或达到最大迭代次数。

四、EM算法的PYhton实现

import numpy as np
from scipy.stats import multivariate_normal

def gaussian_mixture_em(X, n_components, n_iterations=100, tol=1e-4):
    """
    使用EM算法拟合高斯混合模型参数

    参数:
    - X: 输入数据，每行为一个样本
    - n_components: 高斯分布的数量
    - n_iterations: 最大迭代次数
    - tol: 收敛阈值

    返回:
    - means: 各高斯分布的均值
    - covariances: 各高斯分布的协方差矩阵
    - weights: 各高斯分布的权重
    """
    n_samples, n_features = X.shape

    # 初始化参数
    means = X[np.random.choice(n_samples, n_components, replace=False)]
    covariances = [np.eye(n_features) for _ in range(n_components)]
    weights = np.ones(n_components) / n_components

    for _ in range(n_iterations):
        # E步骤
        posterior_probs = np.array([weights[k] * multivariate_normal.pdf(X, means[k], covariances[k]) for k in range(n_components)]).T
        posterior_probs /= posterior_probs.sum(axis=1, keepdims=True)

        # M步骤
        weights = posterior_probs.mean(axis=0)
        means = np.dot(posterior_probs.T, X) / posterior_probs.sum(axis=0, keepdims=True).T
        covariances = [np.dot((X - means[k]).T, (X - means[k]) * posterior_probs[:, k, None]) / posterior_probs[:, k].sum() for k in range(n_components)]

        # 计算似然函数
        log_likelihood = np.sum(np.log(np.sum([weights[k] * multivariate_normal.pdf(X, means[k], covariances[k]) for k in range(n_components)], axis=0)))

        # 判断是否收敛
        if _ > 0 and np.abs(log_likelihood - prev_log_likelihood) < tol:
            break

        prev_log_likelihood = log_likelihood

    return means, covariances, weights

def generate_sample_data():
    """
    生成模拟数据
    """
    np.random.seed(42)
    data1 = np.random.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], 100)
    data2 = np.random.multivariate_normal([5, 5], [[1, -0.8], [-0.8, 1]], 100)
    data = np.vstack([data1, data2])
    return data

def main():
    """
    主函数，用于运行EM算法并输出结果
    """
    # 生成模拟数据
    data = generate_sample_data()

    # 运行EM算法
    n_components = 2
    means, covariances, weights = gaussian_mixture_em(data, n_components)

    # 输出结果
    print("各高斯分布的均值:")
    print(means)
    print("\n各高斯分布的协方差矩阵:")
    print(covariances)
    print("\n各高斯分布的权重:")
    print(weights)

if __name__ == "__main__":
    main()

总结

EM算法是一种迭代优化算法，用于解决含有隐变量的概率模型参数估计问题。其基本思想是通过交替进行期望（Expectation）和最大化（Maximization）两个步骤来优化模型参数。在E步骤中，估计隐变量的条件分布，计算完全数据对数似然的期望值；在M步骤中，通过最大化期望值来更新模型参数。这两步骤交替迭代直至收敛。EM算法广泛应用于聚类、混合模型、隐马尔可夫模型等领域，尤其在处理缺失数据或不完全数据的情况下表现优异。其灵活性和鲁棒性使其成为统计学和机器学习中重要的工具，能够有效解决实际问题中的复杂参数估计挑战。