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{\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反的概率就可以表示为
对上面的两个似然概率进行归一化可以得出后验概率，两者分别是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算法
概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，我们讨论极大似然估计，极大后验概率估计与其类似。
输入：观测变量数据Y，隐变量数据Z，联合分布\(P(Y,Z|\theta)\)，条件分布 \(P(Z|Y,\theta)\)；输出：模型参数\(\theta\)
- (1)选择参数的初值\(\theta^0\)，开始迭代
- (2)E步：记\(\theta^i\)为第\(i\)次迭代参数\(\theta\)的估计值，在第\(i+1\)次迭代的E步，计算
这里，\(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(\theta^{i+1},\theta^{i})-Q(\theta^{i},\theta^{i})|| <\varepsilon_2 \]收敛迭代就结束了。
三、在高斯混合分布中的应用
EM算法的一个重要应用场景就是高斯混合模型的参数估计。高斯混合模型就是由多个高斯模型组合在一起的混合模型（可以理解为多个高斯分布函数的线性组合，理论上高斯混合模型是可以拟合任意类型的分布），例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的：
| 图3 | 图4 |
|---|---|
 |
 |
如果有多个高斯模型，公式表示为：
\[P(y|\theta)=\sum_{k=1}^{K}a_k\phi(y|\theta_{k}) \\ \phi(y|\theta_{k})=\frac{1}{\sqrt{2\pi}\delta_{k}}exp(-\frac{(y-\mu_{k})^2}{2 \delta_{k}^{2}}) \\ a_k>0,\quad \sum a_k =1 \]\(\phi(y|\theta_{k})\)表示为第\(k\)个高斯分布密度模型，定义如上，其中 \(a_k\)表示被选中的概率。在本次模型\(P(y|\theta)\)中，观测数据是已知的，而观测数据具体来自哪个模型是未知的，有点像之前提过的三硬币模型，我们来对比一下，A硬币就像是概率\(a_k\)，用来表明具体的模型，而B、C硬币就是具体的模型，只不过这里有很多个模型，不仅仅是B、C这两个模型。我们用\(\gamma_{jk}\)来表示，则：
\[\gamma_{jk}= \begin{cases} 1& \text{第j个观测数据来源于第k个模型}\\ 0& \text{否则} \end{cases} \]所以一个观测数据\(y_j\)的隐藏数据\((\gamma_{j1},\gamma_{j2},...,\gamma_{jk})\)，那么完全似然函数就是：
\[P(y,\gamma|\theta)= \prod_{k=1}^{K}\prod_{j=1}^{N}[a_{k}\phi(y|\theta_{k})]^{\gamma_{jk}} \]取对数之后等于：
\[\begin{aligned} log(P(y,\gamma|\theta))&=log( \prod_{k=1}^{K}\prod_{j=1}^{N}[a_{k}\phi(y|\theta_{k})]^{\gamma_{jk}})\\ &=\sum_{K}^{k=1}\bigg(\sum_{j=1}^{N}(\gamma_{jk}) log(a_k)+\sum_{j=1}^{N}( \gamma_{jk})\bigg [log(\frac{1}{\sqrt{2\pi}})-log(\delta_{k})-\frac{(y_i-\mu_{k})^2}{2 \delta_{k}^{2}}\bigg]\bigg) \end{aligned} \]- E步 ：
其中我们定义$ \hat{\gamma_{jk}}$：
\[\hat{\gamma\_{jk}} = E(\gamma_{jk}|y,\theta)=\frac{a_k\phi(y_i|\theta_{k})}{\sum_{k=1}^{K}a_k\phi(y_i|\theta_{k}) }\\ j=1,2,..,N；k=1,2,...,K\\ n_k=\sum_{j=i}^{N}E\gamma_{jk} \]于是化简得到：
\[\begin{aligned} Q(\theta，\theta^i) \\ &= \sum_{K}^{k=1}\bigg(n_k log(a_k)+\sum_{j=1}^{N}(E\gamma_{jk})\bigg [log(\frac{1}{\sqrt{2\pi}})-log(\delta_{k})-\frac{(y_i-\mu_{k})^2}{2 \delta_{k}^{2}}\bigg ]\bigg) \end{aligned} \]E步在代码设计上只有\(\hat{\gamma_{jk}}\)有用，用于M步的计算。
- M步 :
对\(Q(\theta,\theta^i)\)求导，得到每个未知量的偏导，使其偏导等于0，求解得到：
\[\hat{m_k}=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}}y_i}{\sum_{j=1}^{N}\hat{\gamma_{jk}}} \\ \hat{\delta_k}=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}}(y_i-\mu_k)^2}{\sum_{j=1}^{N}\hat{\gamma_{jk}}} \hat{a_k}=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}}}{N} \]给一个初始值，来回迭代就可以求得值内容。这一块主要用到了\(Q(\theta，\theta^i)\)的导数，并且用到了E步的\(\hat{\gamma_{jk}}\)。
四、EM算法的PYhton实现
import numpy as np
import random
import math
import time
'''
数据集：伪造数据集（两个高斯分布混合）
数据集长度：1000
------------------------------
运行结果：
----------------------------
the Parameters set is:
alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0
----------------------------
the Parameters predict is:
alpha0:0.4, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9
----------------------------
'''
def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):
'''
初始化数据集
这里通过服从高斯分布的随机函数来伪造数据集
:param mu0: 高斯0的均值
:param sigma0: 高斯0的方差
:param mu1: 高斯1的均值
:param sigma1: 高斯1的方差
:param alpha0: 高斯0的系数
:param alpha1: 高斯1的系数
:return: 混合了两个高斯分布的数据
'''
# 定义数据集长度为1000
length = 1000
# 初始化第一个高斯分布，生成数据，数据长度为length * alpha系数，以此来
# 满足alpha的作用
data0 = np.random.normal(mu0, sigma0, int(length * alpha0))
# 第二个高斯分布的数据
data1 = np.random.normal(mu1, sigma1, int(length * alpha1))
# 初始化总数据集
# 两个高斯分布的数据混合后会放在该数据集中返回
dataSet = []
# 将第一个数据集的内容添加进去
dataSet.extend(data0)
# 添加第二个数据集的数据
dataSet.extend(data1)
# 对总的数据集进行打乱（其实不打乱也没事，只不过打乱一下直观上让人感觉已经混合了
# 读者可以将下面这句话屏蔽以后看看效果是否有差别）
random.shuffle(dataSet)
#返回伪造好的数据集
return dataSet
# 高斯分布公式
def calcGauss(dataSetArr, mu, sigmod):
'''
根据高斯密度函数计算值
依据：“9.3.1 高斯混合模型” 式9.25
注：在公式中y是一个实数，但是在EM算法中(见算法9.2的E步)，需要对每个j
都求一次yjk，在本实例中有1000个可观测数据，因此需要计算1000次。考虑到
在E步时进行1000次高斯计算，程序上比较不简洁，因此这里的y是向量，在numpy
的exp中如果exp内部值为向量，则对向量中每个值进行exp，输出仍是向量的形式。
所以使用向量的形式1次计算即可将所有计算结果得出，程序上较为简洁
:param dataSetArr: 可观测数据集
:param mu: 均值
:param sigmod: 方差
:return: 整个可观测数据集的高斯分布密度（向量形式）
'''
# 计算过程就是依据式9.25写的，没有别的花样
result = (1 / (math.sqrt(2 * math.pi) * sigmod ** 2)) * np.exp(
-1 * (dataSetArr - mu) * (dataSetArr - mu) / (2 * sigmod ** 2))
# 返回结果
return result
def E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):
'''
EM算法中的E步
依据当前模型参数，计算分模型k对观数据y的响应度
:param dataSetArr: 可观测数据y
:param alpha0: 高斯模型0的系数
:param mu0: 高斯模型0的均值
:param sigmod0: 高斯模型0的方差
:param alpha1: 高斯模型1的系数
:param mu1: 高斯模型1的均值
:param sigmod1: 高斯模型1的方差
:return: 两个模型各自的响应度
'''
# 计算y0的响应度
# 先计算模型0的响应度的分子
gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)
# print("gamma0=",gamma0.shape) # 1000, 维向量
# 模型1响应度的分子
gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)
# 两者相加为E步中的分布
sum = gamma0 + gamma1
# 各自相除，得到两个模型的响应度
gamma0 = gamma0 / sum
gamma1 = gamma1 / sum
# 返回两个模型响应度
return gamma0, gamma1
def M_step(muo, mu1, gamma0, gamma1, dataSetArr):
# 依据算法9.2计算各个值
# 这里没什么花样，对照书本公式看看这里就好了
# np.dot 点积：[1,2] [2,3] = [2,6]
mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)
mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)
# math.sqrt 平方根
sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo) ** 2) / np.sum(gamma0))
sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1) ** 2) / np.sum(gamma1))
alpha0_new = np.sum(gamma0) / len(gamma0)
alpha1_new = np.sum(gamma1) / len(gamma1)
# 将更新的值返回
return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new
## 训练主函数
def EM_Train(dataSetList, iter=500):
'''
根据EM算法进行参数估计
算法依据“9.3.2 高斯混合模型参数估计的EM算法” 算法9.2
:param dataSetList:数据集（可观测数据）
:param iter: 迭代次数
:return: 估计的参数
'''
# 将可观测数据y转换为数组形式，主要是为了方便后续运算
dataSetArr = np.array(dataSetList)
# 步骤1：对参数取初值，开始迭代
alpha0 = 0.5
mu0 = 0
sigmod0 = 1
alpha1 = 0.5
mu1 = 1
sigmod1 = 1
# 开始迭代
step = 0
while (step < iter):
# 每次进入一次迭代后迭代次数加1
step += 1
# 步骤2：E步：依据当前模型参数，计算分模型k对观测数据y的响应度
gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)
# 步骤3：M步
mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)
# 迭代结束后将更新后的各参数返回
return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1
if __name__ == '__main__':
start = time.time()
# 设置两个高斯模型进行混合，这里是初始化两个模型各自的参数
# 见“9.3 EM算法在高斯混合模型学习中的应用”
# alpha是“9.3.1 高斯混合模型” 定义9.2中的系数α
# mu0是均值μ
# sigmod是方差σ
# 在设置上两个alpha的和必须为1，其他没有什么具体要求，符合高斯定义就可以
alpha0 = 0.3 # 系数α
mu0 = -2 # 均值μ
sigmod0 = 0.5 # 方差σ
alpha1 = 0.7 # 系数α
mu1 = 0.5 # 均值μ
sigmod1 = 1 # 方差σ
# 初始化数据集
dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)
# 打印设置的参数
print('---------------------------')
print('the Parameters set is:')
print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
))
# 开始EM算法，进行参数估计
alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)
# 打印参数预测结果
print('----------------------------')
print('the Parameters predict is:')
print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (
alpha0, alpha1, mu0, mu1, sigmod0, sigmod1
))
# 打印时间
print('----------------------------')
print('time span:', time.time() - start)
总结
EM算法是一种迭代优化算法，用于解决含有隐变量的概率模型参数估计问题。其基本思想是通过交替进行期望（Expectation）和最大化（Maximization）两个步骤来优化模型参数。在E步骤中，估计隐变量的条件分布，计算完全数据对数似然的期望值；在M步骤中，通过最大化期望值来更新模型参数。这两步骤交替迭代直至收敛。EM算法广泛应用于聚类、混合模型、隐马尔可夫模型等领域，尤其在处理缺失数据或不完全数据的情况下表现优异。其灵活性和鲁棒性使其成为统计学和机器学习中重要的工具，能够有效解决实际问题中的复杂参数估计挑战。
