期望值最大化算法 EM：睹始知终
- 算法思想
- 算法推导
- 算法流程
- E步骤：期望
- M步骤：最大化
- 陷入局部最优的原因
- 算法应用
- 高斯混合模型（Gaussian Mixture Model, GMM）
- 问题描述
- 输入输出
- Python代码实现
算法思想
期望值最大化方法，是宇宙演变、物种进化背后的动力。
如果一个公司在制定年终奖标准时，把每个员工一半的奖金和公司价值观挂钩，人们就会背诵创始人每个语录 — 整个公司都会自动迭代寻找最优解，每个人说话都是公司价值观。
如果一个国家足球不行，把每个孩子的高考分数和足球水平挂钩，人们就会大力投资足球设施，大爷大妈也会把广场让出去给孙子踢足球，谁跟我孙子抢我真的会发疯 — 整个国家都会自动迭代寻找最优解，每个人说话都是公司价值观。
这个思想在算法中就是期望最大化 EM 算法，只要给出一个收益函数， 计算机就会自动的寻找收益最大的那个点。
- 在每一时刻，算出能够最大化收益（期望值）的方向，沿着这个方向走一小步
- 然后再从新的起点重复这个过程，不论从何处起始，最后一定能够达到收益最大的那个终点
EM 算法本质是迭代策略，用于含有隐变量的统计模型中，交替计算期望步骤和最大化步骤，来寻找参数的最优估计。
比如看故事书，但故事中有一些缺失的部分（这些就是隐变量）。
你的目标是填补这些缺失部分，使得整个故事变得连贯和合理。
EM 算法就像一个两步循环过程，帮助你逐渐完善这个故事：
-
期望步骤 (E 步骤)： 在这一步，你根据目前所知的信息，对故事中缺失的部分做出最佳猜测。就好比你根据故事的上下文来推测这些缺失部分可能的内容。
-
最大化步骤 (M 步骤)： 接下来，你根据这些猜测来重新讲述整个故事，并调整故事中其他已知部分的细节，使得整体故事更加合理。这个过程就像根据新的假设来优化故事的连贯性。（M步骤可以使用 MLE 或 MAP）。
这个循环反复进行：你根据当前的故事版本来改善你对缺失部分的猜测，然后再用这些新猜测来优化整个故事。
随着每次迭代，故事变得越来越连贯，直到最终达到一个点，你觉得再怎么调整也无法使故事更好了。
这时，你就找到了最合适的版本来填补缺失部分，你找到了模型参数的最优估计。
再如 市场营销策略：
-
公司在设计营销策略时，通常会试图理解消费者的隐藏需求和偏好（隐藏变量），并据此调整其产品或服务（参数）。
-
通过市场反馈，公司不断调整其策略以最大化销售或品牌影响力，这类似于EM算法的期望步骤和最大化步骤的迭代过程。
算法推导
EM 算法论文：https://web.mit.edu/6.435/www/Dempster77.pdf
概率图模型再复杂都可以简化成俩个变量：观测变量x、隐变量z
比如你正在看一部电影：
- 电影中你能直接看到的场景和角色对话等就像是“观测变量”，这些是你直接获得的信息，不需要猜测或推理。
- 然而，电影也有许多你看不到的部分，比如角色的内心想法、未展示的背景故事，或者导演留下的悬念。这些就像是“隐变量”，你无法直接观察它们，但它们对整个故事的剧情发展（趋势就是人心所向）至关重要。
p ( x ∣ θ ) = ∏ i = 1 n p ( x i ∣ θ ) L ( θ ) = log p ( x ∣ θ ) = ∑ i = 1 n log p ( x i ∣ θ ) = ∑ i = 1 n log ∑ z p ( x i , z i ∣ θ ) \begin{aligned} &p(\mathbf{x}|\theta) \begin{aligned}=\prod_{i=1}^np(x_i|\theta)\end{aligned} \\ &{ L ( \theta )} =\operatorname{log}p(\mathbf{x}|\theta) \\ &=\sum_{i=1}^n\log p(x_i|\theta) \\ &=\sum_{i=1}^n\log\sum_zp(x_i,z_i|\theta) \end{aligned} p(x∣θ)=i=1∏np(xi∣θ)L(θ)=logp(x∣θ)=i=1∑nlogp(xi∣θ)=i=1∑nlogz∑p(xi,zi∣θ)
那我们逐步拆解公式原意：
-
联合概率分布：
- 第一行公式表示，观测数据集 x 在给定参数 θ 的条件下的联合概率分布
- 比如你有 3 张卡片，每张卡片上都有一个秘密数字，这个数字可以是 1、2、3 中的任何一个，我们现在要猜每张卡片上的数字是什么。每张卡片上数字的猜测都是独立的，不会影响其他卡片上的猜测。
- 在数学中，这就是我们说的“联合概率分布”，即我们想知道，所有卡片上每一种可能的数字组合出现的整体概率是多少。
- 如所有卡片上都是1的概率是多少（111）、如所有卡片上是123的概率是多少（123）、（222）、（321）、…、（333） 所有可能的数字组合及其相应的概率。
-
对数似然函数：
- 第二行公式，为了不忘记我们的猜测，我们决定把每次猜的结果写在一个日记本上。因为数字可能很大，所以我们用一种特别的数学“捷径”来记日记，这种捷径就是对数。这样，即使我们猜的数字很大，日记本上的数字也不会太长，更容易计算。
- 在数学中，写在日记本上的这种方法叫做“对数似然函数”，一个帮助我们处理大数字的数学工具。
-
对数似然的求和：
- 第三行公式，现在我们决定把日记本上所有的数字加起来，因为我们用了对数，所以加起来很容易。这就像是玩一个加法游戏，把所有的小数字加起来，得到一个总分。
-
边缘概率：
- 第四行公式，第1张是1、第2张是2，第 3 张卡片藏在盒子里（只有第 3 张未知），我们只知道盒子里可能藏着什么数字（1、2、3）。那先专注于部分已知信息，而忽略未知部分的具体细节，猜对所有看得见的卡片的概率是多少。
- 就是计算 第1张是1、第2张是2 的概率，忽略第三张卡片可能的值。
- 这就是数学中的“边缘概率” —— 它允许我们在部分信息未知的情况下，仍对已知部分进行概率计算。
在概率分布上，就是先猜一个 z 的分布（记为 q），使用 E、M 步骤，去逼近真实分布 L ( θ ) L(\theta) L(θ)：
最后让猜的分布像爬楼梯一样，找到真实分布 L ( θ ) L(\theta) L(θ) 的最高点（最优解）。
用数学公式描述这个过程：
L ( θ ) = ∑ i = 1 n log ∑ z p ( x i , z i ∣ θ ) = ∑ i = 1 n log ∑ z ∞ q i ( z ) p ( x i , z i ∣ θ ) q i ( z ) ≥ ∑ i = 1 n ∑ z q i ( z ) log p ( x i , z i ∣ θ ) q i ( z ) \begin{aligned} L(\theta)& \begin{aligned}=&\sum_{i=1}^n\log\sum_zp(x_i,z_i|\theta)\end{aligned} \\ &\begin{aligned}=&\sum_{i=1}^n\log\sum_z^\infty q_i(z)\frac{p(x_i,z_i|\theta)}{q_i(z)}\end{aligned} \\ &\geq\sum_{i=1}^n\sum_zq_i(z)\log\frac{p(x_i,z_i|\theta)}{q_i(z)} \\ \end{aligned} L(θ)=i=1∑nlogz∑p(xi,zi∣θ)=i=1∑nlogz∑∞qi(z)qi(z)p(xi,zi∣θ)≥i=1∑nz∑qi(z)logqi(z)p(xi,zi∣θ)
-
第一行： L ( θ ) = ∑ i = 1 n log ∑ z p ( x i , z i ∣ θ ) L(\theta) = \sum_{i=1}^n \log \sum_z p(x_i, z_i|\theta) L(θ)=∑i=1nlog∑zp(xi,zi∣θ)
- 比如你正在玩一个寻宝游戏，你有一张地图（ θ \theta θ），地图上标记了很多可能藏宝的地方（这里的藏宝地方就是 x i x_i xi 和 z i z_i zi）。
- x i x_i xi 是你可以看到的地方，而 z i z_i zi 是地图上标记的，但实际上可能藏宝也可能没藏宝的秘密地方。这一行的意思是，你在尝试弄清楚，根据地图，每个地方藏宝的可能性有多大。
-
第二行： = ∑ i = 1 n log ∑ z ∞ q i ( z ) p ( x i , z i ∣ θ ) q i ( z ) = \sum_{i=1}^n \log \sum_z^\infty q_i(z) \frac{p(x_i, z_i|\theta)}{q_i(z)} =∑i=1nlog∑z∞qi(z)qi(z)p(xi,zi∣θ)
- 这一步就像你在用一种特别的放大镜 q i ( z ) q_i(z) qi(z) 来看地图（ θ \theta θ）。
- 这个放大镜可以告诉你，每个秘密地方真的藏宝的机会有多大。
- 你用这个放大镜和地图一起，来计算每个地方可能藏宝的几率。
-
第三行： ≥ ∑ i = 1 n ∑ z q i ( z ) log p ( x i , z i ∣ θ ) q i ( z ) \geq \sum_{i=1}^n \sum_z q_i(z) \log \frac{p(x_i, z_i|\theta)}{q_i(z)} ≥∑i=1n∑zqi(z)logqi(z)p(xi,zi∣θ)
- 最后，这一步就像你在记录你的发现。
- 对于地图上的每一个地方，你都写下了：根据我的放大镜和地图，我认为这里藏宝的机会有多大。”
- 这样，你就得到了一个完整的藏宝地图，上面标记了所有可能藏宝的地方和它们的可能性。
然后根据 Jeasen 不等式，得到公式的下界。
最终的公式是： J ( z , q ) J(z,q) J(z,q)。
- 不断的改变 z，就能不断搜索 θ \theta θ 最大值（概率分布图中的最高点）
于是，EM 算法可分为 E 步骤、M 步骤。
算法流程
E步骤：期望
E 步骤：猜的分布 q 不变，最大化 z。
在图中，q 沿着 x 轴上升，碰到真实分布z 就停止，开始 M 步骤。
M步骤：最大化
M 步骤：猜的分布 q 寻优，z 不变。
在图中，q 沿着 y 轴水平移动，碰不到真实分布z 就停止，开始 E 步骤。
陷入局部最优的原因
EM 算法可能会陷入局部最优。
-
非凸目标函数：EM算法通常用于优化非凸（non-convex）的目标函数。在非凸函数中，可能存在多个局部最优解，这意味着算法可能会在达到一个局部最优点后停止，而这个点不一定是全局最优的。
-
初始值依赖性：EM算法的结果往往依赖于初始参数的选择。如果初始参数选得不好，算法可能会被引导到一个局部最优解而不是全局最优解。
-
迭代方式：EM算法通过交替执行其两个步骤（E步和M步）来逐渐改进参数估计。这种迭代方式可能会导致算法“陷入”某个局部区域的最优解，特别是在目标函数有多个峰值的情况下。
-
模型复杂性和数据的局限性：在一些复杂模型或者数据不足的情况下，EM算法可能无法准确估计出全局最优参数，从而陷入局部最优。
解决这些问题的一种方法是通过多次运行算法，每次使用不同的初始参数，然后从中选择最好的结果。
此外，还可以使用全局优化技术，如模拟退火或遗传算法，来辅助找到更接近全局最优的解。
算法应用
高斯混合模型（Gaussian Mixture Model, GMM）
问题描述
假设我们有一组观测数据点，我们认为这些数据点是由两个不同的高斯分布生成的，但我们不知道每个数据点来自哪个高斯分布。
我们的目标是估计这两个高斯分布的参数（均值和方差）以及每个分布对应的混合系数。
输入输出
- 输入：一组观测数据点。
- 输出：两个高斯分布的参数（均值和方差）和混合系数。
Python代码实现
import numpy as np
from sklearn.mixture import GaussianMixture# 模拟数据生成
np.random.seed(0)
data = np.concatenate([np.random.normal(0, 1, 300), np.random.normal(5, 1.5, 700)]).reshape(-1, 1)# 应用EM算法
gmm = GaussianMixture(n_components=2, random_state=0)
gmm.fit(data)# 输出结果
print(f'均值: {gmm.means_.ravel()}')
print(f'方差: {gmm.covariances_.ravel()}')
print(f'混合系数: {gmm.weights_.ravel()}')
这段代码首先生成了一些模拟数据，数据是由两个不同的高斯分布合成的。
然后使用sklearn库中的GaussianMixture模型来应用EM算法。最后，打印出两个高斯分布的均值、方差以及混合系数。
