1977年，DempSter首次提出EM算法。
假设四种实验结果，发生的概率依次为 ，且发生的次数为 ，求 的估计。
解：使用MLE，得到：
上式是关于 的一元三次方程，不易解。
因此，以下另作处理（引入隐变量）：
将第一部分 分为 ，且出现次数为 次
将第三部分 分为 ，且出现次数为 次；
则
(1)
现在，并不知道 （隐变量）的值，只能知道分布的信息， 服从的分布为二项分布，概率数值类似于条件概率，第一个的概率是用 除以得到的，第二个同理：
其中， ，
第一步（E步）：求期望的目的是为了消去隐变量 。
代入(1)式，得到：
第二步（M步）：取最大值。
EM算法使用迭代法来更新参数。 （精髓）
任意取 ，就可以开始按照上面的公式进行迭代了。
收敛性 ：
DempSter证明：在很一般的条件下，最后会收敛。（可以参考李航老师的《统计学习方法》）
解析解：能列出公式解决的，数值上是更准确的（相比迭代解），比如MLE就是列出公式求解。
迭代解：退而求其次，当解析解难求的时候，通过迭代逼近的方式，可以获得令人满意的解，比如EM就是为了解决当MLE遇到高次方程难以求解的时候，提出的方法。
问：给定参数 ，观测变量 ，隐变量 ，如何估计参数 ？
从观测序列，可以获得：
此时，对数似然函数为：
由于包含和（积分）的对数，因此直接求解困难。
解析解困难，转而使用迭代解：假设第i次迭代后的为 ，由于我们希望似然函数 是增大的，即 。
此时，考虑两者的差：
不等式右边是 的下界，记为 ，那么，使得下界尽可能大，即：
Algorithm: Estimation Maximum (EM)
举例：以三硬币模型为例。有A、B、C三枚硬币，分别有 的概率为正面。每次试验为：先投A硬币，如果A为正面，则投B硬币；否则，投C硬币。最终，可以观测到的结果为硬币的正/反面，但是不知道是由B还是C投出的（隐变量）。问：如果某次试验数为10的结果为：{1,1,0,1,0,0,1,0,1,1}，如何估计参数 ？
显然，题目的 隐变量为A硬币投出的结果，此时可以采用EM解法。
先从“E”入手，求解Q函数：
然后，逐一击破：
回代 函数：
极大似然求导数，令其为0，能取得极值点：
令上式为0
------对应书(9.6)式
令上式为0
------对应书(9.7)式
令上式为0
------对应书(9.8)式
至此，只要根据当前迭代下的 ，就能得到不同下标的 ，进而得到下一次迭代的 。
r语言中bpinom函数的基本用法为：pbinom(x,size,,prob)，该函数为事件的累积概率，它用于表示概率的单个值。例如：抛掷硬币100次，正面向上不超过50次的概率，即pbinom(50,100,0.5)。
r语言有四个内置函数来生成二项分布。它们的描述分别如下：
dbinom(x,size,prob)函数，该函数表示每个点的概率密度分布。
pbinom(x,size,prob)函数，该函数为事件的累积概率，它表示概率的单个值。
qbinom(p,size,prob)函数，该函数采用概率值，并给出累积值与概率值匹配的数字。
rbinom(n,size,prob)函数，该函数从给定样本产生给定概率的所需数量的随机值。
其中，x是数字的向量，p是概率向量，n是观察的数量，size是试验的数量，prob是每个试验成功的概率。
最近在做文本挖掘的时候遇到了EM算法，虽然读书的时候简单地接触过，但当时并没有深入地去了解，导致现在只记得算法的名字。既然EM算法被列为数据挖掘的十大算法之一，正好借这个机会，重新学习一下这个经典的算法。学习的过程中，我发现网上的资料大多讲解地不够细致，很多地方解释得并不明了。因此我决定抛开别人的想法，仅从数学推导本身出发，尽力理解每一个公式的含义，并将其对应到实际的实验过程当中。这篇博客记录了我对与EM算法的思考与理解，也是我人生中的第一篇博客，希望能够对于想要学习EM算法的同学有所帮助。
前面谈到我在做文本挖掘的时候遇到了EM算法，EM算法用于估计模型中的参数。提到参数估计，最常见的方法莫过于极大似然估计——在所有的候选参数中，我们选择的参数应该让样本出现的概率最大。相信看到这篇笔记的同学一定对极大似然估计非常熟悉，而EM算法可以看作是极大似然估计的一个扩充，这里就让我们用极大似然估计来解决一个简单的例子，来开始正式的讨论。
有A，B，C三枚硬币，我们想要估计A，B，C三枚硬币抛出正面的概率 , , 。我们按如下流程进行实验100次：
记录100次实验的结果如下：
我们将上面的实验结果表述如下：
表示第i次实验中，硬币A的结果，1代表正面，0代表反面； 表示第i次实验中，硬币B或硬币C抛出正面的个数，则参数 的极大似然估计分别为：
即硬币A，B，C各自抛出正面的次数占总次数的比例，其中 为指示函数。
实验流程与1相同，但是我们不慎遗失了硬币A的记录结果，导致我们只知道随后十次抛出了多少次正面，多少次反面，却不知道实验结果来自于硬币B还是硬币C。在这种情况下，我们是否还能估计出 , , 的值呢？
这时候利用极大似然估计似乎行不通了， 因为这种情况下，我们不但缺失了硬币A产生的观测值，同时也不知道哪些观测值属于硬币B，哪些观测值属于硬币C。
有些同学可能会提出，虽然我们无法得到三个硬币各自产生的样本，但是我们依然可以得到每个观测值出现的概率。比如在第一次实验中， 我们抛出了5次正面5次反面，我们可以做如下思考：
假设这5次正面由硬币B得到，那么概率应该为 ，而这次观测值来自于硬币B，也就是硬币A抛出正面的概率为
假设这5次正面由硬币C得到，那么概率应该为 ，而这次观测值来自于硬币C，也就是硬币A抛出反面的概率为
综合起来，利用条件概率公式，这个观测值出现的概率就是
因此我们可以将样本整体的概率和似然函数利用 , , 表示出来，通过对似然函数求导，令其关于 的偏导数等于0，我们可以求出三个参数的值。
这个思路听上去十分合理，我们可以顺着这个思路进行数学推导，看看可以得到什么样的结果。首先我们计算样本的概率:
对应的似然函数为
其中 关于 的条件分布为
的分布为
因此我们可以得到
至此，我们成功地得到了似然函数。然而观察可以发现，这个函数是由100项对数函数相加组成，每个对数函数内部包含一个求和，想通过求导并解出导数的零点几乎是不可能的。当然我们可以通过梯度下降来极小化这个函数，借助深度学习库的自动微分系统在实现上也非常容易。但是这种做法过于简单粗暴，有没有办法来优雅地解决这个问题呢？在继续讨论之前，我们先将这类问题进行一般化表述：
我们观测到随机变量 产生的m个相互独立的样本 , 的分布由联合分布 决定， 是缺失数据或无法在实验中被直接观测到，称为 隐变量 ，我们想要从样本中估计出模型参数 的值。在接下来的讨论中，我们假定 的取值是离散的，于是可以得到似然函数如下:
接下来，我们就探讨一下，如何利用EM算法解决这个问题。
这一部分的数学推导，主要参考了吴恩达CS229n的笔记，并且根据个人的思考和理解，尽力对公式的每一步进行详细的解释。我们先简单地介绍一下琴生不等式。
琴生不等式有多种形式，下面给出其离散形式的表述和概率论中的表述:
1.若 为严格凹函数， 为定义域内的n个点， 是n个正实数，且满足 , 则下述不等式成立:
当且仅当 时，不等式取等号。
2.若 为严格凹函数， 为实值随机变量，且期望存在，则下述不等式成立:
当且仅当 ，即 为常数时，不等式取等号。
注： 这里将函数上方为凹集的函数称为凹函数， 例如 函数就是凹函数。
相信大家对琴生不等式都十分熟悉，因此这里就不做过多的说明。接下来，我们将琴生不等式应用到我们的问题中。
回到我们之前的问题上， 我们想要极大化下面这个函数:
但是我们无法对这个函数直接求导，因此我们借助琴生不等式，对这个函数进行变换。为了让过程看上去简洁，下面只对求和中的第 项进行计算。
令 满足 ，且 ，则根据琴生不等式，可以得到：
当且仅当 为常数时，上述不等式取等号。也就是说，对于任意 ， 是一个与 无关的量。设对于任意 ，我们可以得到：
因此当 时，不等式 取等号，容易验证此时 , 与 无关。将 综合一下，我们可以得到以下结论:
到这里为止，我们已经拥有了推导出EM算法的全部数学基础，基于 我们可以构建出E步和M步。上面的数学推导虽然看上去略为复杂，但实际上只用到了三个知识点：
1.琴生不等式:
2.条件概率:
3.联合分布求和等于边缘分布:
对上面的数学推导有疑问的同学，可以结合上面这三点，再将整个推导过程耐心地看一遍。
大部分关于EM算法的资料，只是在数学形式上引入了 函数，即 ，以满足琴生不等式的使用条件，却没有过多地解释 函数本身。这导致了很多人完全看懂了算法的推导，却还是不理解这些数学公式究竟在做什么，甚至不明白EM算法为什么叫做EM算法。所以在给出E步和M步之前，我想先谈一谈 函数。
我们回顾一下 函数所满足的条件（暂时不考虑琴生不等式取等号的限制），
在 所有可能的取值处有定义。可以看出， 是 的样本空间上任意的一个概率分布。因此，我们可以对不等式 进行改写。首先我们可以将含有 的求和写成期望的形式:
这里 指的是在概率分布 下，求随机变量 和 的期望。有同学会问，为什么我们平时求期望的时候只要写 ，并没有指明是在哪个概率分布下的期望。这是因为一般情况下，我们都清楚地知道随机变量 所服从的分布 ，并且默认在分布 下求期望。
举个例子，我手上有一个硬币，抛了10次，问抛出正面次数的期望。这种情况下，大部分人会默认硬币是均匀的，也就是说抛出正面的次数 服从二项分布 ，期望 。这时有人提出了质疑，他说我认为你这个硬币有问题，抛出正面的概率只有0.3，那么在他眼里， 期望 。
回到正题，我们利用等式 改写不等式 ，可以得到:
这正是琴生不等式在概率论中的形式。我们可以将不等式倒过来理解：
首先，假定随机变量 服从概率分布 ， 是 的样本空间上的任意一个概率分布。这里 可以是一组定值，也可以是关于参数 的函数。
显然，当我们取不同的 时，随机变量 的期望也会随之改变。需要注意的是，由于 与 相关，所以这里的期望不是一个数值，而是关于 的函数。
当我们令 为 的后验分布 时，上面的期望最大。这里有两点需要注意，1. 后验分布 也是一个关于参数 的函数。2. 由于期望是关于 的函数，所以这里的最大指的并非是最大值，而是最大的函数。
若对于每一个 ，我们都令 为 的后验分布 ，则上述期望之和等于我们要极大化的似然函数，即
通过上述分析，我们为寻找似然函数的极大值点 提供了一个思路。我们不去极大化似然函数本身，而是去极大化 。至于如何将这个思路实际应用，就要利用到EM算法中的E-step和M-step。
这一节中，我们先给出E-step和M-step的数学形式，随后在结合抛硬币的例子来解释这两步究竟在做什么。下面进入算法的流程，首先我们任意初始化 ，按下述过程进行迭代直至收敛：
在第 次迭代中，
(E-step)对于每个 ，令
(M-step)更新 的估计值，令
EM算法从任意一点 出发，依次利用E-step优化 ，M-step优化 ，重复上述过程从而逐渐逼近极大值点。而这个过程究竟是怎样的呢，就让我们一步步地揭开EM算法的面纱。
假设我们现在随机初始化了 ，进入第一轮迭代：
(E-step)
由于我们已经假定模型参数为 ，所以此时 不再是与 有关的函数，而是由一组常数构成的概率分布。结合抛硬币的例子来看，这一步是在我们已知模型参数 的基础上(虽然这是我们瞎猜的)，去推测每一次的观测值是由哪个硬币产生的，或者说我们对每一次观测值做一个软分类。比如我们根据初始化的参数，计算出 ， 。可以解释为第 个观测值有20%的概率来自于硬币B，80%的概率来自于硬币C；或者说硬币A抛出了0.2个正面，0.8个反面。
(M-step)
考虑到 是一组常数，我们可以舍弃常数项，进一步简化上面这个要极大化的函数
由于 不再与 相关，因此上面的函数变成了对数函数求和的形式，这个函数通常来说是容易求导的，令导数等于0，我们可以求出新的参数 。我们仍旧以抛硬币为例进行解释，
令 , 可以得到，
这三个参数的解释是显而易见的。我们在E-step中对每个观测值进行了软分类， 可以看成是硬币A抛出正面的次数，所以 是 的极大似然估计； 是我们抛硬币B的次数， 是硬币B抛出正面的次数，所以 是 的极大似然估计；对于 我们有相同的解释。
我们将这个结果与抛硬币1中极大似然估计的结果相比较可以发现，之前结果中的指示函数 变成了这里的 ，在指示函数下，某个观测值要么来自于硬币B，要么来自于硬币C，因此也称为硬分类。而在 函数下，某个观测值可以一部分来自于硬币B，一部分来自于硬币C，因此也称作软分类。
将上述两步综合起来，EM算法可以总结如下：我们首先初始化模型的参数，我们基于这个参数对每一个隐变量进行分类，此时相当于我们观测到了隐变量。有了隐变量的观测值之后，原来含有隐变量的模型变成了不含隐变量的模型，因此我们可以直接使用极大似然估计来更新模型的参数，再基于新的参数开始新一轮的迭代，直到参数收敛。接来下我们就讨论为什么参数一定会收敛。
前面写了太多的公式，但是这一部分我不打算给出收敛性的数学推导。其实数学上证明EM算法的收敛性很容易，只需要证明每一轮迭代之后，参数的似然函数递增，即
