使用最大期望(EM)算法与贝叶斯定理解决GMM问题

Expectation-Maximization Algorithm (EM) 算法与贝叶斯定理

• 注1: 本文仅基于GMM高斯混合模型的EM算法解法来讨论

• 注2: 本文不会涉及证明EM算法对于期望似然函数 (ELF) 的 weak monotonic increase 性质，但这个证明是这个算法可行的基础，请自行证明

---

问题引入:

• 我们现在有两个已知服从正态分布的数据集混在一起了（比如最经典的，两个班级的成绩数据混在一起了），我们要给出一个合理的预测每个成绩更可能是那个班的
• 我们首先创建两个基于正太分布的数据集，我们现在假设班级1有60个人，班级2有40个人（这个数据后面有用）。
• 接下来我会生成一组数据，并在数轴上展示:
  • 班级 1 的真实参数： $num_1, \mu_1, \sigma_1 = 50, 80, 5$
  • 班级 2 的真实参数： $num_2, \mu_2, \sigma_2 = 40, 50, 10$

---

分析问题:

高斯混合模型GMM

• 首先，这是一个高斯混合模型(GMM),即多个高斯分布混合在一起
• 我们的目的是找到最合适的两个高斯分布来分解这个分布
• 即: 优化两个高斯分布的 四个参数 使得预测分布与真实分布的的期望似然函数最大（这个期望似然函数可以简单理解为loss function）
• 注意，EM算法的根本目的是找到四个参数使得期望似然函数最大

为什么要用EM：

• EM算法用于在存在隐藏或未观察到的变量的情况下，估计概率模型参数。
• EM算法在存在隐变量的情况下相对于其他优化算法，比如梯度下算法，更有优势。因为我们不需要知道“隐变量”具体是什么

在这个情景中与EM的关系:

• 依据前文， EM算法需要有两个要素：隐变量 和 参数， 而在这个模型中:
  • 隐变量: 这个数据是一班的还是二班的？我们可以通过设隐变量 Z = 1 或 2 表示
  • 参数: 两个正态分布，共需要四个变量来表示

---

EM 算法简述:

初始化:

我们将随机或按一定规律初始化四个参数

E步骤

• 1. 给出当前似然函数，即估计的四个参数表示的两个正态分布的pdf表达式
• 2. 根据当前的似然函数求出 后验概率

M步骤

• 根据后验概率更新计算新的四个参数
• 重新进入E步骤开始下一步迭代

正确性证明: 略

---

从机器学习与统计而不是概率的视角看待贝叶斯定理

既然涉及到计算 “后验概率”， 那么我们就绕不开贝叶斯定律

在更一般的解释框架下，后验概率（Posterior Probability） 是表示在观察到某些数据 $x$ 之后，某个假设或模型 $\theta$ 是真实的概率。与似然函数一样，后验概率也用于量化一个假设或模型与观测数据的一致性，但它还进一步地结合了先验信息。

用数学语言表示，后验概率 $P(\theta | x)$ 通常通过贝叶斯定理计算：

$$P(\theta | x) = \frac{P(x | \theta) \cdot P(\theta)}{P(x)}$$

其中：

• $P(\theta)$ 是先验概率，表示在观察数据之前我们对 $\theta$ 的信念或知识。

• $P(x | \theta)$ 是似然函数， 表示在给定 $\theta$ 的条件下观测到数据 $x$ 的概率或相对可能性。 似然函数的意义如下：

  • 在贝叶斯公式中，似然概率 $P(x | \theta)$ 量化了给定模型或假设 $\theta$ 的情况下，观察到数据 $x$ 的相对可能性或概率。。 高似然值意味着在给定模型或假设 $\theta$ 的情况下，数据 $x$ 出现的概率更高，这通常会导致这个值在该模型或假设的后验概率中也更高。
  • 总体而言，似然在贝叶斯分析中提供了一种量化观测数据与模型之间一致性或拟合程度的方法。通过与先验概率 $P(\theta)$ 和边缘概率 $P(x)$ 的结合，似然帮助我们更新和修正对模型或假设 $\theta$ 的信念，从而达到根据数据进行推理的目的。

• $P(x)$ 是边缘概率（也称为“证据”），

  • 表示 在考虑所有可能的模型或假设 $\theta$ 后，观察到数据 $x$ 的总概率 $P(x)$
  • 反映了数据 $x$ 在没有任何先验信息或假设限制的情况下的“自然”出现概率。
  • 由于evidence是所有可能的模型和假设概率和，所以也当作归一化分母，对应的，分子是只考虑了一种模型或者假设下观察到某个数据的概率
  • 这有助于我们理解在给定数据 $x$ 的情况下，不同的模型或假设 $\theta$ 的相对可能性。
  • 通常通过对所有可能的 $\theta$ 值进行积分或求和来计算：$P(x) = \int P(x | \theta) \cdot P(\theta) \, d\theta$（连续情况）或 $P(x) = \sum_{\theta} P(x | \theta) \cdot P(\theta)$（离散情况）。

• 后验概率 $P(\theta | x)$

  • 综合了先验概率和观测数据，为我们提供了一个更新的、更全面的视图
  • 用于评估假设或模型 $\theta$ 的可信度。

贝叶斯在模型预估中的作用

贝叶斯定理计算的后验概率提供了一种统一的方式，用于结合先验信息和新的观测数据，以更新我们对某个假设或模型的信念。

定义RV与事件:

• 我们定义 Random Variable $X$ 与 Random Variable $Z$
  • $X$ 表示这个样本的成绩，依此定义事件 $X = x_i \ ( x_i ∈ [0, 100] )$
  • $Z$ 表示这个样本是哪个班级的，依此定义事件 $Z = z_i \ ( z_i ∈ \set{1, 2} )$
  • 注意， $X$, $Z$ 两个是一个样本的两个“成员变量”，关系是等价的

带入贝叶斯定理

• 先验概率 prior: $p(Z = z_i)$ 表示如果没有任何其他条件的话一个样本属于某个班级的概率是什么
• 似然 likelihood: $p(X = x_i | Z = z_i)$
  • 这个似然表示的是 如果这个样本属于z_i班的话，这个成绩出现的概率， 就是z_i班的人拿这个成绩的概率
  • 但是实际上这个似然有什么意义呢？ 这个似然其实是表示了现在的模型对这个成绩归属的估计到底靠不靠谱？
  • 人话：如果一班平均分非常非常高，那么对于一个具有非常低的成绩的样本，我们理应认为他不太可能是一班的，但是如果先验概率认为他是一班的，那么就是先验概率的问题
• 证据（边缘概率）evidence: $p(X =x_i)$ 表示在没有任何其他条件的情况下，这个样本得到这个成绩的概率
  • 计算方法应该是: $\sum_{i = 1}^n{p(Z = x_i | Z = z_i)p(Z = z_i)}$
• 后验概率 posterior: $p(Z = z_i | X = x_i)$ 表示如果这个样本的成绩是x_i的话，这个样本是z_i班级的可能性
• 公式:

$$p(Z = z_i | X = x_i) = \frac{p(X = x_i | Z = z_i)p(Z = z_i)}{\sum_{i = 1}^n{p(X = x_i | Z = z_i)p(Z = z_i)}}$$

---

EM 算法步骤:

知道了贝叶斯定理之后就可以围绕贝叶斯定理来设计出基于贝叶斯定理的EM算法

---

E:

• 首先对于E，要求是计算出后验概率
• 我们就把所有上文的数据带入以下公式计算出所以p(Z = z_i | X = x_i):

$$p(Z = z_i | X = x_i) = \frac{p(X = x_i | Z = z_i)p(Z = z_i)}{\sum_{i = 1}^n{p(X = x_i | Z = z_i)p(Z = z_i)}}$$

---

M

• M步骤就是想办法基于E中算出的后验概率来更新似然函数的参数来最大化期望似然函数

首先， 如何更新μ才能比原先的μ更接近真实的μ呢？（μ值的intuition对于σ值也是类似的）

• intuition():
  • 我们要构造一个函数，使得这个函数可以“慢慢把$\mu$值向正确的方向拖动”，实际上我们会构造以下函数（本文未给出严禁证明）:

$$\mu_{new}(1) = \sum_{i = 1}^n{\frac{p(Z = 1 | X = x_i) \cdot x_i}{\sum_{i = 1}^n{p(Z = 1 | X = x_i)}}}$$

或者如果把第i个数据的后验概率定义为$\omega_i$可以简写为:

$$\mu_{new}(1) = \sum_{i = 1}^n{\frac{\omega_i \cdot x_i}{\sum_{i = 1}^n{\omega_i}}}$$

理解: 这个公式每个component作用如下:
分子：$\sum_{i = 1}^n{\omega_i \cdot x_i}$, 含义为 $x_i$对于后验概率的加权和（注意不是加权平均）
分母: $\sum{i = 1}^n{\omega_i}$, 含义为归一化，让加权和变成加权平均！！

为什么会越来越接近呢（定性分析）

• 首先，明确一点，如果数据量足够大，那么真实中心附近一定会有越来约密集的样本点
• 如果$\mu$偏离了真实值，那么，原本真实中心的$\omega$后验概率权重会偏小（这是似然函数偏小导致的），然后估计的中心样本点比较稀疏，但是给了非常大的权重。
• 这时候算法的加权平均会慢慢把$\mu$往真实值那里拉，如以下例子:
  • 我给估计的样本中心的1个元素权重定为5，实际样本中心的10个元素的权重定为1，最后加权平均一定会偏向实际样本中心
  • 但是这样不可保证一定会到“最大似然函数”，但是能非常好的接近

不严谨证明这个方法的可行性:

首先我们现在想要证明的命题是:

$$当\mu恰好等于\mu_{true}时，这个公式计算出来的\mu_{new}刚好就是\mu_{true}$$

• 证明:

$$如果数据样本足够大，则\mu_{true} = \sum_{i = 1}^n{p(X = x_i)x_i}$$

• 假设 $p(Z = 1 | X = x_i) = p(X = x_i)$

$$\mu_{new}(1) = \sum_{i = 1}^n{\frac{\omega_i \cdot x_i}{\sum_{i = 1}^n{\omega_i}}} = \sum_{i = 1}^n{\frac{p(X = x_i) \cdot x_i}{\sum_{i = 1}^n{p(X = x_i)}}} = \sum_{i = 1}^n{p(X = x_i) \cdot x_i} = \mu_{true} \ (补充：\sum_{i = 1}^n{p(X = x_i)} = 1)$$

• 所以仅需证明，在n足够大时:

$$p(Z = 1 | X = x_i) = p(X = x_i)$$

• 证明：

$$p(X = x_i) = p(X = x_i | Z = 1)p(Z = 1) + p(X = x_i | Z = 2)p(Z = 2)$$

• 其中由于分布就是真实分布，所以$p(Z = 1), p(z = 2)$中有一个为0一共为1，假设$Z = 1$, 有：

$$p(X = x_i) = p(X = x_i | Z = 1)p(Z = 1)$$

• 根据贝叶斯定理

$$p(Z = z_i | X = x_i) = \frac{p(X = x_i | Z = z_i)p(Z = z_i)}{\sum_{i = 1}^n{p(X = x_i | Z = z_i)p(Z = z_i)}}$$

• 由于就是真实分布，所以Evidence和为1，故有:

$$p(Z = 1 | X = x_i) = p(X = x_i | Z = 1)p(Z = 1) = p(X = x_i)$$

Q.E.D

这个公式对于σ²的扩展:

与$\mu$类似，方差也有类似的迭代公式，只是加权平均的不再是样本的成绩, 而是 $ (x_i - \mu_{new})^2 $ , 样本成绩与新的中心的距离的平方差，公式如下:

$$\sigma^2_{new}(1) = \sum_{i = 1}^n{\frac{\omega_i \cdot (x_i - \mu_{new})^2}{\sum_{i = 1}^n{\omega_i}}}$$

注: 这里$\omega_i$表示第i个样本的后验概率，与上文同

先验概率的迭代公式（已知占比）：

与$\mu, \sigma^2$不同的是，先验概率prior的初始值是确定的，而不是随机的
prior的公式是:$p(Z = z_i)$, 表明，没有任何条件的情况下对于任意的样本应该如何估计他的归属
根据定义可以非常容易得知，这个东西就是一个班级的人数占比啊！
所以计算公式理所应当的为:

$$prior_1 = \frac{num_1}{num_1 + num_2}$$

$$prior_2 = \frac{num_2}{num_1 + num_2}$$

注意，如果已知每个班的人数，这个可以当作一个参数不变

先验概率迭代公式（未知占比）：

这样初值就只能用随机了
intuition: 既然不知到占比，那么就假设我们的后验概率是准确的，后验概率的意义是$p(Z = z_i | X = x_i)$，即每个点属于某个班级的概率， 那么就直接去后验概率的平均，表示平均每个点分到某个班的概率是多少就可以了，公式也会非常简单:

$$prior_1 = \frac{1}{n}\sum_{i = 1}^n{\omega_i}$$

注: 这里$\omega_i$表示第i个样本的后验概率，与上文同

---

实现核心代码

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.stats import norm

# 初始化参数
init_mu1, init_sigma1 = 85, 8  # 班级 1
init_mu2, init_sigma2 = 55, 9  # 班级 2
init_pi1 = num1 / (num1 + num2)  # 最开始的先验概率priority，即班级 1 的学生所占比例
init_pi2 = 1 - init_pi1     # 班级 2 的学生所占比例

# 参数存储
mu1, sigma1 = init_mu1, init_sigma1
mu2, sigma2 = init_mu2, init_sigma2
pi1, pi2 = init_pi1, init_pi2
iterations = 10

# 存储每次迭代的参数以便后续可视化
mu1_list, sigma1_list = [mu1], [sigma1]
mu2_list, sigma2_list = [mu2], [sigma2]

for iter in range(iterations):
    # E步: 计算后验概率（贝叶斯公式）
    # 首先计算先验概率，即最开始设定的p(X = x_i | Z = 1)
    likelihood1 = norm.pdf(data, mu1, sigma1)
    likelihood2 = norm.pdf(data, mu2, sigma2)
    
    # 计算每个数据点属于每个班级的未标准化概率，即p(X = xi | Z = 1)p(Z = 1)
    numerator1 = pi1 * likelihood1
    numerator2 = pi2 * likelihood2
    
    # 计算归一化因子（证据）
    evidence = numerator1 + numerator2
    
    # 计算标准化后的后验概率
    posterior1 = numerator1 / evidence
    posterior2 = numerator2 / evidence
    
    # M步: 更新参数
    mu1 = np.sum(posterior1 * data) / np.sum(posterior1)
    sigma1 = np.sqrt(np.sum(posterior1 * (data - mu1)**2) / np.sum(posterior1))
    pi1 = np.sum(posterior1) / len(data)
    
    mu2 = np.sum(posterior2 * data) / np.sum(posterior2)
    sigma2 = np.sqrt(np.sum(posterior2 * (data - mu2)**2) / np.sum(posterior2))
    pi2 = np.sum(posterior2) / len(data)
    
    # 存储参数
    mu1_list.append(mu1)
    sigma1_list.append(sigma1)
    mu2_list.append(mu2)
    sigma2_list.append(sigma2)

---

结果可视化: