EM(Expectation-maximization algorithm)译为期望最大化算法，EM算法是数据挖掘的十大算法之一，主要解决有隐含变量时，如何利用最大似然法求解未知参数。现实中会遇到多个类别数据混杂在一起，各个类别数据虽然是一个概率分布，但数学期望或方差不同，每次取得一个数据时并不知道是归于哪个类别下，样本数据属于哪个类别的信息是一个隐含变量，遇到这种情况时不能直接用最大似然法。

    一、EM算法解决抛硬币问题

    1.1 了解EM算法过程

    根据大数定律，扔硬币时只要次数足够多，硬币的正反的概率应该各为50%,或者换个说法，真币正面的概率应该是50%。现有两个假币A和B,两种硬币的朝正面的概率未知（肯定不是50%），测试人员每次测试时从两枚假币中任意选一枚扔10次为一组数据，共得到5组数据，记录硬币正反的情况，由于实验中，实验每一组数据时没做标记，所以实验者不知道每组是硬币A还是硬币B，要得到A和B朝正面的概率，首先需要知道每组实验中硬币是A或B的最大可能性分布，此分布对于利用最大似然法求解而言是隐含变量。

    再来看另外一个例子：两类正态分布的数据混合在一起，其方差一样，数学期望不一样，现在需要分别求两类数据的数学期望。在用最大似然法求解参数时，在这个场景下每个数据都有两种可能性，数据的可能性也可以理解为一个分布，即这组混合数据中两类数据各占百分之多少，这是需要先求解出来的隐含变量。

    回到抛硬币的问题上来，这个问题是一个二项分布问题，EM算法解决这类问题方法分两步：

E步骤：要求解未知参数θA,θB，先随机初始化未知参数θA,θB，利用θA,θB求出每一组数据的分布,在硬币问题中，该分布指每一组实验中硬币分别是A和B的概率。

M步骤：求出每组实验数据分布后，综合每组数据求出一套新的参数θA,θB，与初始参数比较，如果发现误差在设定范围内则停止，否则拿新的参数继续E步骤。

    可以看出，E步骤目的是求解隐含变量从而获得实验数据分布，得到分布后，隐变量此时是一个常数，这时就可以使用最大似然法了，M步骤实质是利用最大似然法求目标参数值。EM算法解决硬币问题可以用下图来表示:

    上图a代表每次知道扔的是硬币A或者是B，这种场景下由于没有隐含变量，可以用传统最大似然法求出A和B的'朝正’概率，譬如硬币A一共被扔了30次，其中24次是正面，得到硬币A正面的概率是80%，类似可以得到硬币B正面的概率是45%，这两个参数可以验证接下来EM算法得到概率是否准确。

    图b代表每次并不知道扔的是哪个硬币，此时用的EM算法，算法首先初始化一个概率，假设A的正面概率为0.6,B的正面概率是0.5，这两个概率事件是独立不相关的，所以其和并不等于1。接下来选择其中一个硬币扔10次为一组数据；拿第一组为例，扔出的结果是[H,T,T,T,H,H,T,H,T,H],正面5次反面5次，假设是A货币，由θA=0.6并结合二项概率分布知，出现这组情况的概率是：

0.65*(1-0.6)5=0.0007962624

而如果是硬币B,可以计算得到出现这组数据的概率是：

0.55*(1-0.5)5=0.0009765625

把上面分布归一化处理得到第一组数据分别为A和B的概率：

P(A)=0.0007962624/(0.0007962624＋0.0009765625)≈0.45

P(B)=0.0009765625/(0.0007962624＋0.0009765625)≈0.55

在先设定θA=0.6，θB=0.5后出现[H,T,T,T,H,H,T,H,T,H]的情况时，是A币的概率是0.45，是B的概率是0.55；根据这个分布还可以得到分别为A或B时的正反面期望次数,例如第一组数据中共出现了5次正面，5次反面:

A币正面的期望次数：5(5次正面)*0.45=2.2次 

B币正面的期望次数：5(5次正面)*0.55=2.8次 

5次为反的情况可以这样计算:

A币反面的期望次数：5(5次反面)*0.45=2.2次 

B币反面的期望次数：5(5次反面)*0.55=2.8次   

每组数据是A币或是B币是一个互斥事件，概率之和是1，这里需要和θA、θB区别一下。这样就得到第一组数据分布情况，按此方法可同样计算其他四组数据分布，以上的过程是EM算法的E步骤。

    在得到每组数据分布之后，接下来是M步骤：分别用A、B正面'期望次数'除以正反面总次数，得到新的θA、θB, 判断新的参数与原参数的误差是否在可接受范围内，如果相差较大，则利用新的参数重复之前操作；如果误差比较小则结束计算，从而得到最终解。E和M步骤一前一后像一套组合拳，在EM算法中先E后M为一次迭代过程。

1.2  推导EM算法过程

EM算法就是围绕公式(1.4)展开，最大似然法是求能使得公式(1.4)获得最大值时的参数θ，但这里多出来一组隐含变量zj，上面提到，zj是一个随机变量：条件概率pθ(xi|zj)代表当设定硬币为A或B后，出现该组数据的概率。在示例中，设定硬币A正面概率为0.6后，第一组数据有pθ(x1|z1)=0.65*(1-0.6)5=0.0007962624，同样对于硬币B有pθ(x1|z2)=0.55*(1-0.5)5=0.0009765625；pθ(zj)代表每组数据分别是硬币A和B的概率分布，示例中是通过将pθ(xi|z1)和pθ(xi|z2)归一化后求得的。通过示例可以看出，求解pθ(zj)、pθ(xi|zj)依赖参数θ，而欲求θ必须知道A和B的分布即每轮zi值，这便形成了一个求解的死循环，可利用Jensen不等式将隐含随机变量zi变为常数:

上图是y=ln(x)的函数图像，用x表示CD之间任意一个变量，E是CD中点可表示为数学期望E(x)，E点在函数图像对应F点，函数值为f(E(x))；点H是割线AB的中点，它的纵坐标是CD两端函数值的平均，所以点H纵坐标可用E(f(x))表示,观察图形有:

f(E(x))≥E(f(x))    (1.5)

需要关注(1.5)式何时取等号，通过观察发现当AB这段曲线是直线时，F、H将完全重合，此时(1.5)式取等号，有了这个思路后，引入一个新函数:

     (1.6)

从形式上看Qi(zj)归一化后和值为1，Qi(zj)含义为第i组数据中硬币为A或B的概率分布，将Qi(zj)代入到(1.4)式中并结合Jensen不等式(1.5)有以下关系式:

上面推导过程使用了换元法，可将理解为公式(1.5)的x,(1.5)式中f(x)换元后变为，相应的(1.5)中的f(E(x))换元后变为，这里用新生成函数Qi(zj)代表概率分布，类似的有：

上面推导公式可得到下面的关系式:

 （1.7）

通过之前分析，(1.7)式取等号时，必须为线性函数，对于对数函数而言，只有是常数才满足此条件，设在t轮迭代初始参数为θt,比如示例中第一轮迭代有初始参数θ1=(θA,θB)=(0.6,0.5)：

第一组数据为硬币A时概率有

第一组数据为硬币B时概率有

这说明有了初始参数θt后，pθt(xi,zj)是一个常数值，pθt(xi,zj)代表设定初始参数后，硬币分别为A、B时可得到该组数据结果的概率，有了这个先决条件后，可设，C为常数：

上式推导中最后一步利用了贝叶斯概率公式简化了形式，至此得到一个结论：当有关系式时，公式(1.7)等式成立，而Qi(zj)含义是第i组数据中硬币为A或B的概率分布，在示例中是通过归一化求出的：

综上所述，求解公式（1.4）最优解时，由于存在隐含变量无法使用最大似然法，而当公式(1.7)取等号时，可以将隐含变量变为常数，此时就可以用最大似然法求未知参数θ值了。也可以从函数角度来分析这个过程，利用Jensen不等式得到目标函数的下界函数：

为了能求函数最大值，就需要通过提升下界函数来实现，这是水涨船高的原理。当公式(1.7)取等号时代表下界函数与相切于一点，此时已到最高位置，不能再上升了，再上升便不再是下界函数了，而公式(1.7)取等号成立时，通过之前分析可以知道，必须为一个常数，从这个条件又可以推导出隐含变量Qi(zj)值，将隐含变量变为常量、消掉隐藏变量后，变为一个仅有参数θ函数：

 （1.8）

此时对求导，利用最大似然法就可以得到本轮的参数θt+1，求导取最大值过程对应于EM算法的M步骤，以上分析过程可以通过下图表示：

上面红色线代表目标函数，紫色代表下界函数。下界函数通过E步骤逐步提升接近红色线，蓝色线代表下界函数与目标函数相切于一点，即公式(1.7)取相等时的情形；接下来通过M步骤：求蓝色曲线最大值所对应θ，得到本轮迭代的参数θt ＋1，这时判断蓝色曲线最大值与目标函数即红色线最大值的误差，如果误差在接受范围内结束迭代，否则用新的θt ＋1参数开始下一轮EM算法。如果把上面这个三维图像做一个横切面即可得到一个二维图像：

1.2 利用Python代码实现EM算法

针对以上问题给出一个python代码示例