EM算法
EM算法，全称Expectation Maximization Algorithm，译作最大期望化算法或期望最大算法，它是一种迭代算法，用于含有隐变量（hidden variable）的概率参数模型的最大似然估计或极大后验概率估计。
结论：
E-step：计算期望
\[\mathcal{Q}\big(\theta;\theta^{(t)}\big)=\mathbb{E}_{p(\mathbf{z}\mid \mathbf{x};\mathbf{\theta}^{(t)}}\big[\log p(\mathbf{x},\mathbf{z};\mathbf{\theta})\big] =\int_Z P(Z\mid X,\theta^{(t)}) \log P(X,Z\mid \theta) dZ\\ PPT中以离散为例子，所以是\sum_{z}P(Z\mid X,\theta^{(t)})替换\int_Z P(Z\mid X,\theta^{(t)})\]M-step:更新参数
\[\boldsymbol{\theta}^{(t+1)}=\arg\max_\theta Q(\boldsymbol{\theta};\boldsymbol{\theta}^{(t)})\]总述
期望最大（Expectation Maximization）算法是一种从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。EM算法是机器学习十大算法之一，或许确实是因它在实际中的效果很好吧。
定义
EM算法，全称Expectation Maximization Algorithm，译作最大期望化算法或期望最大算法，它是一种迭代算法，用于含有隐变量（hidden variable）的概率参数模型的最大似然估计或极大后验概率估计。
例子
感性例子A：
例子简介：
假设现在有两枚硬币1和2，,随机抛掷后正面朝上概率分别为$P_1$，$P_2$。为了估计这两个概率，做实验，每次取一枚硬币，连掷5下，记录下结果，如下：
| 硬币 | 结果 | 统计 |
|---|---|---|
| 1 | 正正反正反 | 3正-2反 |
| 2 | 反反正正反 | 2正-3反 |
| 1 | 正反反反反 | 1正-4反 |
| 2 | 正反反正正 | 3正-2反 |
| 1 | 反正正反反 | 2正-3反 |
可以很容易地估计出$ P_1$和$ P_2$，（硬币1实验15次正面6次，硬币2实验10次正面5次）如下：
\[P_1 = ( 3 + 1 + 2 ) / 15 = 0.4 \\ P_2 = ( 2 + 3 ) / 10 = 0.5\]到这里，一切似乎很美好，下面我们加大难度。
加入隐变量z
还是上面的问题，现在我们抹去每轮投掷时使用的硬币标记，如下：
| 硬币 | 结果 | 统计 |
|---|---|---|
| Unknown | 正正反正反 | 3正-2反 |
| Unknown | 反反正正反 | 2正-3反 |
| Unknown | 正反反反反 | 1正-4反 |
| Unknown | 正反反正正 | 3正-2反 |
| Unknown | 反正正反反 | 2正-3反 |
好了，现在我们的目标没变，还是估计$ P_1$和$ P_2$，要怎么做呢？
显然，此时我们多了一个隐变量$z$，可以把它认为是一个5维的向量$(z_1,z_2,z_3,z_4,z_5)$，代表每次投掷时所使用的硬币，比如$ z_1$，就代表第一轮投掷时使用的硬币是1还是2。但是，这个变量$ z$不知道，就无法去估计$P_1$和$P_2$，所以，我们必须先估计出$ z$，然后才能进一步估计$P_1$和$P_2$。
但要估计$z$，我们又得知道$ P_1$和$ P_2$，这样我们才能用最大似然概率法则去估计$ z$，这不是鸡生蛋和蛋生鸡的问题吗，如何破？
答案就是先随机初始化一个$ P_1$和$P_2$，用它来估计$ z$，然后基于$ z$，还是按照最大似然概率法则去估计新的$P_1$和$ P_2$，如果新的$P_1$和$P_2$和我们初始化的$ P_1$和$ P_2$2一样，请问这说明了什么？
这说明我们初始化的$ P_1$和$ P_2$是一个相当靠谱的估计！
就是说，我们初始化的$ P_1$和$ P_2$，按照最大似然概率就可以估计出$z$，然后基于$z$，按照最大似然概率可以反过来估计出$ P_1$和$ P_2$，当与我们初始化的$ P_1$和$ P_2$一样时，说明是$ P_1$和$ P_2$很有可能就是真实的值。这里面包含了两个交互的最大似然估计。
如果新估计出来的$ P_1$和$ P_2$和我们初始化的值差别很大，怎么办呢？就是继续用新的$ P_1$和$ P_2$迭代，直至收敛。
这就是下面的EM初级版。
EM初级版
我们不妨这样，先随便给$ P_1$和$ P_2$赋一个值，比如：
\[P_1 = 0.2\\ P_2 = 0.7\]然后，我们看看第一轮抛掷最可能是哪个硬币。 如果是硬币1，得出正正反正反的概率为（注意这里是有先后顺序的）
\[0.2*0.2*0.2*0.8*0.8 = 0.00512\]如果是硬币2，得出正正反正反的概率为
\[0.7*0.7*0.7*0.3*0.3=0.03087\]然后依次求出其他4轮中的相应概率。做成表格如下：
| 轮数 | 若是硬币1 | 若是硬币2 |
|---|---|---|
| 1 | 0.00512 | 0.03087 |
| 2 | 0.02048 | 0.01323 |
| 3 | 0.08192 | 0.00567 |
| 4 | 0.00512 | 0.03087 |
| 5 | 0.02048 | 0.01323 |
按照最大似然法则： 第1轮中最有可能的是硬币2 第2轮中最有可能的是硬币1 第3轮中最有可能的是硬币1 第4轮中最有可能的是硬币2 第5轮中最有可能的是硬币1
我们就把上面的值作为$z$的估计值。然后按照最大似然概率法则来估计新的$ P_1$和$ P_2$。
我们就把上面的值作为z zz的估计值。然后按照最大似然概率法则来估计新的$P_1$和$P_2$。
\[P_1 = ( 2 + 1 + 2 ) / 15 = 0.33 \\ P_2 = ( 3 + 3 ) / 10 = 0.6\]设想我们是全知的神，知道每轮抛掷时的硬币就是如前面例子中的那样，$P_1$和$P_2$的最大似然估计就是0.4和0.5（下文中将这两个值称为$P_1$和$P_2$的真实值）。那么对比下我们初始化的$P_1$和$P_2$和新估计出的$P_1$和$P_2$：
| 初始化的P1 | 估计出的P1 | 真实的P1 | 初始化的P2 | 估计出的P2 | 真实的P2 |
|---|---|---|---|---|---|
| 0.2 | 0.33 | 0.4 | 0.7 | 0.6 | 0.5 |
看到没？我们估计的$ P_1$和$ P_2$相比于它们的初始值，更接近它们的真实值了！
可以期待，我们继续按照上面的思路，用估计出的$ P_1$和$ P_2$再来估计$z$，再用$z$来估计新的$ P_1$和$ P_2$，反复迭代下去，就可以最终得到$P_1 = 0.4，P_2=0.5$，此时无论怎样迭代，$ P_1$和$ P_2$的值都会保持0.4和0.5不变，于是乎，我们就找到了$ P_1$和$ P_2$的最大似然估计。
这里有两个问题： 1、新估计出的$P_1$和$P_2$一定会更接近真实的$P_1$和$P_2$？ 答案是：没错，一定会更接近真实的$P_1$和$P_2$，数学可以证明，但这超出了本文的主题，请参阅其他书籍或文章。
2、迭代一定会收敛到真实的$ P_1$和$ P_2$吗？ 答案是：不一定，取决于$P_1$和$P_2$的初始化值，上面我们之所以能收敛到$P_1$和$P_2$，是因为我们幸运地找到了好的初始化值。
EM进阶版
下面，我们思考下，上面的方法还有没有改进的余地？
我们是用最大似然概率法则估计出的$z$值，然后再用$z$值按照最大似然概率法则估计新的$ P_1$和$ P_2$。也就是说，我们使用了一个最可能的$z$值，而不是所有可能的$z$值。
如果考虑所有可能的$z$值，对每一个$z$值都估计出一个新的$ P_1$和$ P_2$，将每一个$z$值概率大小作为权重，将所有新的$ P_1$和$ P_2$分别加权相加，这样的$ P_1$和$ P_2$应该会更好一些。
所有的$z$值有多少个呢？显然，有$2^5=32$种，需要我们进行32次估值？
不需要，我们可以用期望来简化运算。
| 轮数 | 若是硬币1 | 若是硬币2 |
|---|---|---|
| 1 | 0.00512 | 0.03087 |
| 2 | 0.02048 | 0.01323 |
| 3 | 0.08192 | 0.00567 |
| 4 | 0.00512 | 0.03087 |
| 5 | 0.02048 | 0.01323 |
利用上面这个表，我们可以算出每轮抛掷中使用硬币1或者使用硬币2的概率。比如第1轮，使用硬币1的概率是：
\[0.00512 / ( 0.00512 + 0.03087 ) = 0.14\]使用硬币2的概率是
\(1 − 0.14 = 0.86\) 依次可以算出其他4轮的概率，如下：
| 轮数 | $ z_i$=硬币1 | $z_i$=硬币2 |
|---|---|---|
| 1 | 0.14 | 0.86 |
| 2 | 0.61 | 0.39 |
| 3 | 0.94 | 0.06 |
| 4 | 0.14 | 0.86 |
| 5 | 0.61 | 0.39 |
上表中的右两列表示期望值。看第一行，0.86表示，从期望的角度看，这轮抛掷使用硬币2的概率是0.86。相比于前面的方法，我们按照最大似然概率，直接将第1轮估计为用的硬币2，此时的我们更加谨慎，我们只说，有0.14的概率是硬币1，有0.86的概率是硬币2，不再是非此即彼。这样我们在估计$P_1$或者$P_2$时，就可以用上全部的数据，而不是部分的数据，显然这样会更好一些。
这一步，我们实际上是估计出了$z$的概率分布，这步被称作E步。
结合硬币 A 的概率和下面的投掷结果，我们利用期望可以求出硬币 A 和硬币 B 的贡献
结合下表：
| 硬币 | 结果 | 统计 |
|---|---|---|
| Unknown | 正正反正反 | 3正-2反 |
| Unknown | 反反正正反 | 2正-3反 |
| Unknown | 正反反反反 | 1正-4反 |
| Unknown | 正反反正正 | 3正-2反 |
| Unknown | 反正正反反 | 2正-3反 |
我们按照期望最大似然概率的法则来估计新的$ P_1$和$ P_2$：
以$P_1$估计为例，这里我们要做的事情是在之前推算出来的使用硬币1的概率下，去计算出现现在这种抛掷情况的可能性，就是说计算一下期望看看有多少次等价的硬币1抛出来为正，多少次抛出来为负，以第一轮为例，硬币1在第一轮中被使用的该概率是0.14，现在出现了三次正面两次反面。
第1轮的3正2反相当于
\[0.14∗3=0.42正\\ 0.14∗2=0.28反\]依次算出其他四轮，列表如下：
| 轮数 | 正面 | 反面 |
|---|---|---|
| 1 | 0.42 | 0.28 |
| 2 | 1.22 | 1.83 |
| 3 | 0.94 | 3.76 |
| 4 | 0.42 | 0.28 |
| 5 | 1.22 | 1.83 |
| 总计 | 4.22 | 7.98 |
上面这个表格这样理解：也就是coin A的情况下出现正面的期望和出现反面的期望
| 轮数 | coin A 正面 反面 | coin B 正面 反面 |
|---|---|---|
| 1 | 0.42 0.28 | 略 |
| 2 | 1.22 1.83 | |
| 3 | 0.94 3.76 | |
| 4 | 0.42 0.28 | |
| 5 | 1.22 1.83 | |
| 总计 | 4.22 7.98 |
可以看到，改变了$z$值的估计方法后，新估计出的$P_1$要更加接近0.4。原因就是我们使用了所有抛掷的数据，而不是之前只使用了部分的数据。
这步中，我们根据E步中求出的$z$的概率分布，依据最大似然概率法则去估计$ P_1$和$ P_2$，被称作M步。
例子B
如果说例子 A 需要计算你可能没那么直观，那就举更一个简单的例子：
现在一个班里有 50 个男生和 50 个女生，且男女生分开。我们假定男生的身高服从正态分布： $N(\mu_1,\sigma^2_1)$ ，女生的身高则服从另一个正态分布： $N(\mu_2,\sigma^2_2)$ 。这时候我们可以用极大似然法（MLE），分别通过这 50 个男生和 50 个女生的样本来估计这两个正态分布的参数。
但现在我们让情况复杂一点，就是这 50 个男生和 50 个女生混在一起了。我们拥有 100 个人的身高数据，却不知道这 100 个人每一个是男生还是女生。
这时候情况就有点尴尬，因为通常来说，我们只有知道了精确的男女身高的正态分布参数我们才能知道每一个人更有可能是男生还是女生。但从另一方面去考量，我们只有知道了每个人是男生还是女生才能尽可能准确地估计男女各自身高的正态分布的参数。
这个时候有人就想到我们必须从某一点开始，并用迭代的办法去解决这个问题：我们先设定男生身高和女生身高分布的几个参数（初始值），然后根据这些参数去判断每一个样本（人）是男生还是女生，之后根据标注后的样本再反过来重新估计参数。之后再多次重复这个过程，直至稳定。这个算法也就是 EM 算法
例子总结
以上，我们用一个实际的小例子，来实际演示了EM算法背后的idea，共性存于个性之中，通过这个例子，我们可以对EM算法究竟在干什么有一个深刻感性的认识，掌握EM算法的思想精髓。
推导
以下来自深入理解EM算法（ELBO+KL形式），加入个人理解，另一种理解是用 Jensen不等式，网上大部分是用Jensen不等式来推导的
还有另一个很详细的教程：EM算法与GMM（高斯混合聚类）
期望最大化(Expectation Maximization, EM)算法一般用于解决带有隐藏数据的概率模型的参数估计问题。那么什么是隐藏数据呢？举个栗子，比如现在我们拿到一堆身高数据，这里面既有男生的身高也有女生的身高，但是具体每个数据是来自男生还是女生我们是不知道的，这里每个数据所对应的性别就属于隐藏数据。
给定一组观测数据 $ X = x_i,i=1,…N $ ,以及模型参数 $ \theta$ ，假设对应的隐藏数据 $ Z= {z_i}_{i=1}^N$ ，我们称 $ (x_i,z_i)$ 为完整数据。如果不考虑隐藏数据，我们就可以直接使用极大似然估计的方法估计出参数 $ \theta$ :
\[\theta_{MLE} = \arg\max\limits_\theta\log P(X\mid \theta) = \arg\max\limits_\theta\sum_i^N\log P(x_i\mid \theta)~~(1) \\\]但由于隐藏数据的存在， 我们有
\[P(X\mid \theta) = \int_Z P(X, Z\mid \theta)dZ~~(2) \\\]可以看到，如果将(2)带入(1)中，在log里面会出现积分(求和)符号，导致对似然函数的求导变得困难，无法求解。对于这种无法直接求解的问题，我们通常会采用迭代求解的策略，一步一步逼近最终的结果，在EM算法中就是E步和M步的交替进行，直至收敛。
EM算法
在推导EM算法之前，我们先给出它的算法步骤。模型参数被初始化为 $ \theta^{(0)} $，假设在第 t 次迭代后，参数的估计值为 $ \theta^{(t)} $，对于第 t+1 次迭代，具体分为两步：
- E-step：求期望
- M-step: 最大化 $Q(\theta,\theta^{(t)})$ 并求解$ \theta^{(t+1)}$
推导 (ELBO+KL形式)
tips：下面符号记法$\mid 和;$有混用，一般右边是$\theta$的话就可以认为等价
根据条件概率公式，我们有 $P(X;\theta) = \frac{P(X,Z;\theta)}{P(Z;X, \theta)}$
两边同时取对数，有$ \log P(X;\theta) = \log P(X,Z;\theta) - \log P(Z\mid X,\theta)$
下面的构造就比较有技巧性了。我们引入 $Z $的概率分布为 $q(Z)$ ，并且有$ \int_Zq(Z)dZ = 1$（连续）或$\sum_{z}q(z)=1$（离散）。我们将上式构造为
\[\log P(X;\theta)=\sum_{Z}q(Z)\log\frac{P(X,Z;\theta)q(Z)}{\underbrace{q(Z)P(Z \mid X;\theta)}}\\ =\sum_{Z}q(Z)\log \frac{P(X,Z;\theta)}{q(Z)} + \sum_{Z}q(Z)\log \frac{q(Z)}{P(Z\mid X,\theta)} \\\]理解方法一：
上式中第二项即为KL散度的定义，于是：
\[\log P(X;\theta)=L(q;\theta) + KL(q\mid P(Z\mid X,\theta))\]当q的分布和P相同，也就是$q(Z)=P(Z\mid X,\theta)$时，有$KL(q\mid P(Z\mid X,\theta)=0$，此时$\log P(X;\theta)=L(q;\theta)$，所以
\[\log P(X;\theta)=L(q;\theta)=\sum_{Z}q(Z)\log \frac{P(X,Z;\theta)}{q(Z)}\\ =\sum_{Z}P(Z\mid X,\theta)\log \frac{P(X,Z;\theta)}{P(Z\mid X,\theta)}\\ =\sum_{Z}P(Z\mid X,\theta)\log P(X,Z;\theta)-\sum_{Z}P(Z\mid X,\theta)P(Z\mid X,\theta)\]后者为常数，最大化时可以去掉，所以目标最大化
\[Q(\theta,\theta^{(t)})=\sum_{Z}P(Z\mid X,\theta)\log P(X,Z;\theta)=E_{P(Z\mid x;\theta^{(t)})}(logP(X,Z;\theta))\]理解方法二：
拆开后即为：
\[\log P(X;\theta) = \log \frac{P(X,Z;\theta)}{q(Z)} - \log \frac{P(Z\mid X,\theta)}{q(Z)} \\\]然后两边同时求关于变量 Z 的期望
\[\mathbb{E}_Z[\log P(X;\theta)] = \mathbb{E}_Z[\log \frac{P(X,Z;\theta)}{q(Z)}] - \mathbb{E}_Z[\log \frac{P(Z\mid X,\theta)}{q(Z)}] \\\]将期望写成积分的形式
\[\int_Z q(Z)\log P(X;\theta)dZ = \int_Zq(Z)\log \frac{P(X,Z;\theta)}{q(Z)}dZ - \int_Zq(Z)\log \frac{P(Z\mid X,\theta)}{q(Z)}dZ \\\]同时由于$ \log P(X\mid \theta) $和 $Z$ 无关，所以
\[\log P(X;\theta) = \int_Zq(Z)\log \frac{P(X,Z;\theta)}{q(Z)}dZ - \int_Zq(Z)\log \frac{P(Z\mid X,\theta)}{q(Z)}dZ \\\]注意最右边的积分项 $- \int_Zq(Z)\log \frac{P(Z\mid X,\theta)}{q(Z)}dZ $，这个其实就是$ q(Z)和 P(Z\mid X,\theta) $之间的相对熵，即$ Kullback-Leibler divergence (KL divergence)$，记作
\[D_{KL}(q(Z)\mid \mid P(Z\mid X,\theta)) = \int_Zq(Z)\log \frac{q(Z)}{P(Z\mid X,\theta)}dZ \\\]所以我们有
\[\log P(X;\theta) = \int_Zq(Z)\log \frac{P(X,Z;\theta)}{q(Z)}dZ + D_{KL}(q(Z)\mid \mid P(Z\mid X,\theta)) \\\]根据KL divergence的性质，我们知道 $D_{KL}(q(Z)\mid \mid P(Z\mid X,\theta)) \ge 0$ 当且仅当 $q(Z)=P(Z\mid X,\theta)$取等号，因此
\[\log P(X;\theta) \ge \int_Zq(Z)\log \frac{P(X,Z;\theta)}{q(Z)}dZ \\\]我们便得到了 $\log P(X;\theta)$ 的一个下界 ，我们称之为Evidence Lower Bound (ELBO)，后面就可以通过迭代的方式不断抬高ELBO，使得 $\log P(X;\theta) $增大。但目前还有一个问题，就是 $q(Z)$是未知的，所以下界还是没法求，怎么办呢？我们可以直接在每轮迭代时令 $q(Z)=P(Z\mid X,\theta^{(t)}) $，此时$ D_{KL}(q(Z)\mid \mid P(Z\mid X,\theta^{(t)})) = 0$ ，因为我们想要ELBO和 $\log P(X;\theta) $的差距尽可能小，这样抬高ELBO才会使得 $\log P(X;\theta) $增益更大，所以将KL这一项直接置为0是比较合理的，此时ELBO就变为
\[\int_ZP(Z\mid X,\theta^{(t)})\log \frac{P(X,Z;\theta)}{P(Z\mid X,\theta^{(t)})}dZ = \mathbb{E}_{Z\mid X,\theta^{(t)}}[\log\frac{P(X,Z\mid \theta)}{P(Z\mid X,\theta^{(t)})}] \\\]展开有
\[\mathbb{E}_{Z\mid X,\theta^{(t)}}[\log\frac{P(X,Z\mid \theta)}{P(Z\mid X,\theta^{(t)})}] = \mathbb{E}_{Z\mid X,\theta^{(t)}}[\log P(X,Z\mid \theta)] - \mathbb{E}_{Z\mid X,\theta^{(t)}}[\log P(Z\mid X,\theta^{(t)}] \\\]因为我们最终的目标是求出某个$ \hat\theta $使得ELBO最大，上式的第二项与 $\theta$ 无关，可看成是一个常数，所以可以直接扔掉也不会影响$ \hat\theta $的值，那么式子变成了
\[\mathbb{E}_{Z\mid X,\theta^{(t)}}[\log P(X,Z\mid \theta)] \\\]这样我们得到了最开始所提到的EM算法E-step求期望的那个式子。紧接着就是求解 $\theta $,使得该期望达到最大，即M-step
\[\theta^{(t+1)} = \arg\max\limits_\theta \mathbb{E}_{Z\mid X,\theta^{(t)}}[\log P(X,Z\mid \theta)] \\\]以上便是EM算法的ELBO+KL形式的推导过程，当然还有其他的推导方式，后续我们再做介绍。
总结
EM算法是一种迭代式求解的方法，通过重复E-step和M-step, 逐步抬高证据下界来最大化目标函数的值，直至收敛到某个最优点。需要注意的是，EM算法一般收敛到局部最优，无法得到全局最优解。另外，在机器学习领域中一个最常见的使用到了EM思想的算法便是K-means聚类：每个聚类簇的中心是我们需要估计的参数，每个样本的所属类别可看作是隐藏数据。对于E步，我们根据上一轮的聚类结果对类中心进行更新；对于M步，每个样本又被重新聚类到这些更新后的类中心里面。重复E步和M步，直到类中心不再变化，由此完成了对样本的K-means聚类。
GMM-EM算法步骤总结
机器学习包括无监督学习；无监督学习包括聚类；GMM-EM算法是一个聚类算法，它将EM算法应用于GMM，适用于凸性数据的聚类。
以下步骤建议先去理解GMM算法后再来看，更容易理解，点这里跳转到GMM算法介绍
算法步骤总结如下：
- 初始化参数
- E step 更新后验分布
$\gamma_{nk}$ ，意思是数据 $n$（数据$n$即为第$n$个数据）符合分布$ k$(分布$k$即为第$k$个簇内的分布)的概率
- M step 更新参数
- 检查是否收敛，否则令
随后返回 E step。
其中定义了后验分布（特别地，这里称作“责任（responsibility）”；
\[\gamma_{nk}=p(z_{nk}=1 \mid \pmb x_n)\quad (n=1,\cdots,N;k=1,\cdots,K)\\\]可以参考《聚类》中的GMM算法：
$\gamma_{ji}$ ，意思是数据 $j$（数据j即为第j个数据）符合分布$ i$(分布i即为第i个簇内的分布)的概率
还定义了各类别的“有效数量”
\[N_k = \sum_{n=1}^N \gamma_{nk}\quad (k=1,\cdots,K)\\\]- 使用EM算法求得 GMM 的最大似然解之后，可以更进一步，对每个 $\pmb x_n $，估计出它最可能来自哪个类别
至此完成了凸性数据的无监督聚类。
其中，后验概率更新公式和参数更新公式推导如下。
参考资料
【机器学习笔记10】EM算法——直观理解与详细推导_似然函数等于1-CSDN博客
