硬币~隐藏~鸡蛋~期望
一、EM算法实例讲解
首先本文参考：https://www.jianshu.com/p/1121509ac1dc
一个简单的例子：
假设现在有两枚硬币1和2，随机抛掷后正面朝上概率分别为P1，P2。为了估计这两个概率，做实验，每次取一枚硬币，连掷5下，记录下结果，如下：
可以很容易地估计出P1和P2，如下：(P1 = （3+1+2）/ 15 = 0.4,P2= （2+3）/ 10 = 0.5)，到这里，一切似乎很美好，下面我们加大难度
加入隐变量z：
还是上面的问题，假设我们不知道每轮使用的是哪个硬币，如下：
好了，现在我们的目标没变，还是估计P1和P2，要怎么做呢？
显然，此时我们多了一个隐变量z，可以把它认为是一个5维的向量（z1,z2,z3,z4,z5)，代表每次投掷时所使用的硬币，比如z1，就代表第一轮投掷时使用的硬币是1还是2。但是，这个变量z不知道，就无法去估计P1和P2，所以，我们必须先估计出z，然后才能进一步估计P1和P2。
但要估计z，我们又得知道P1和P2，这样我们才能用最大似然概率法则去估计z，这不是鸡生蛋和蛋生鸡的问题吗，如何破？
答案就是先随机初始化一个P1和P2，用它来估计z，然后基于z，还是按照最大似然概率法则去估计新的P1和P2，如果新的P1和P2和我们初始化的P1和P2一样，请问这说明了什么？（此处思考1分钟）
这说明我们初始化的P1和P2是一个相当靠谱的估计！就是说，我们初始化的P1和P2，按照最大似然概率就可以估计出z，然后基于z，按照最大似然概率可以反过来估计出P1和P2，当与我们初始化的P1和P2一样时，说明是P1和P2很有可能就是真实的值。这里面包含了两个交互的最大似然估计。
如果新估计出来的P1和P2和我们初始化的值差别很大，怎么办呢？就是继续用新的P1和P2迭代，直至收敛。
这就是下面的EM初级版。
EM初级版
我们不妨这样，先随便给P1和P2赋一个值，比如：P1 = 0.2，P2 = 0.7，然后，我们看看第一轮抛掷最可能是哪个硬币。
如果是硬币1，得出3正2反的概率为：0.2*0.2*0.8*0.2*0.8 = 0.00512
如果是硬币2，得出3正2反的概率为：0.7*0.7*0.3*0.7*0.3 = 0.03087
然后依次求出其他4轮中的相应概率。做成表格如下：
按照最大似然法则：
第1轮中最有可能的是硬币2
第2轮中最有可能的是硬币1
第3轮中最有可能的是硬币1
第4轮中最有可能的是硬币2
第5轮中最有可能的是硬币1
我们就把上面的值作为z的估计值。然后按照最大似然概率法则来估计新的P1和P2。
P1 =（2+1+2）/ 15 = 0.33
P2 =（3+3）/ 10 = 0.6
设想我们是全知的神，知道每轮抛掷时的硬币就是如本文第001部分标示的那样，那么，P1和P2的最大似然估计就是0.4和0.5（下文中将这两个值称为P1和P2的真实值）。那么对比下我们初始化的P1和P2和新估计出的P1和P2：
看到没？我们估计的P1和P2相比于它们的初始值，更接近它们的真实值了！
可以期待，我们继续按照上面的思路，用估计出的P1和P2再来估计z，再用z来估计新的P1和P2，反复迭代下去，就可以最终得到P1 = 0.4，P2=0.5，此时无论怎样迭代，P1和P2的值都会保持0.4和0.5不变，于是乎，我们就找到了P1和P2的最大似然估计
我们再试一轮，看看P1和P2是否更接近0.4和0.5，上一轮我们已经知道P1=0.33，P2=0.6，我们再次计算出每一轮的频率：
按照最大似然法则：
第1轮中最有可能的是硬币2
第2轮中最有可能的是硬币1
第3轮中最有可能的是硬币1
第4轮中最有可能的是硬币2
第5轮中最有可能的是硬币1
我们就把上面的值作为z的估计值。然后按照最大似然概率法则来估计新的P1和P2：
P1 =（2+1+2）/ 15 = 0.33
P2 =（3+3）/ 10 = 0.6
可以看出已经收敛了，再进行下一轮，硬币序列不会变，P1和P2的值也不会变
这里有两个问题：
1、新估计出的P1和P2一定会更接近真实的P1和P2？
答案是：没错，一定会更接近真实的P1和P2，数学可以证明，但这超出了本文的主题，请参阅其他书籍或文章。
2、迭代一定会收敛到真实的P1和P2吗？
答案是：不一定，取决于P1和P2的初始化值，上面我们之所以能收敛到P1和P2，是因为我们幸运地找到了好的初始化值。
EM进阶版
下面，我们思考下，上面的方法还有没有改进的余地？
我们是用最大似然概率法则估计出的z值，然后再用z值按照最大似然概率法则估计新的P1和P2。也就是说，我们使用了一个最可能的z值，而不是所有可能的z值。
我们可以用期望来简化运算
利用上面这个表，我们可以算出每轮抛掷中使用硬币1或者使用硬币2的概率。比如第1轮：
使用硬币1的概率是：0.00512 / (0.00512+0.03087) = 0.14
使用硬币2的概率是：1 - 0.14 = 0.86
依次可以算出其他4轮的概率，如下：
上表中的右两列表示期望值。看第一行，0.86表示，从期望的角度看，这轮抛掷使用硬币2的概率是0.86。相比于前面的方法，我们按照最大似然概率，直接将第1轮估计为用的硬币2，此时的我们更加谨慎。
我们只说，有0.14的概率是硬币1，有0.86的概率是硬币2，不再是非此即彼。这样我们在估计P1或者P2时，就可以用上全部的数据，而不是部分的数据，显然这样会更好一些。
这一步，我们实际上是估计出了z的概率分布，这步被称作E步。
结合下表：
我们按照期望最大似然概率的法则来估计新的P1和P2：以P1估计为例，第1轮的3正2反相当于
0.14 * 3 = 0.42正
0.14 * 2 = 0.28反
依次算出其他四轮，列表如下：
P1=4.22 / (4.22 + 7.98) = 0.35
可以看到，改变了z值的估计方法后，新估计出的P1要更加接近0.4。原因就是我们使用了所有抛掷的数据，而不是之前只使用了部分的数据。
这步中，我们根据E步中求出的z的概率分布，依据最大似然概率法则去估计P1和P2，被称作M步。
二、EM算法及推广
我们再来看看<<李航*统计学习>>第九章中关于EM算法介绍：含有隐变量的概率模型参数的极大似然估计法（这个隐含变量就是对应的上面例子中的两个硬币，不知道选择的是哪个硬币，所以称为隐变量）
如：在一次独立实验中，有三枚硬币A，B，C，正面出现的概率分别是：((pi,p,q))进行如下投掷实验：先投掷A，若是正面就选B，若是反面就选C，然后最终结果正面出现记为1，反面记为0。独立重复n次实验。（如n=10），观测结果如下：
((1,1,0,1,0,0,1,0,1,1))
假设我们只能看到最终的观测结果，不能观测到投掷硬币的过程，问如何估计三硬币正面出现的概率，即求((pi,p,q))对应的值
一次独立实验模型可以写作：
(P(y| heta) = sum_zP(y,z| heta) = sum_zP(z| heta)P(y|z, heta) = pi p^y(1-p)^{1-y} + (1-pi)q^y(1-q)^{1-y}) ((9.1))
其中y表示一次观测变量的结果：1或者0。z表示隐变量，表示未观测到的硬币A的结果，z的结果不可观测。( heta=(pi,p,q))是模型参数。
由于进行了n次独立重复实验，每次都会产生一个隐变量结果和观测结果，其中我们将观测数据表示为(Y = (Y_1,Y_2,...,Y_n)^T)，未观测到的数据表示为(Z = (Z_1,Z_2,...,Z_n)^T)，则整个观测过程的似然函数为：
(P(Y|Theta) = sum_ZP(Z|Theta)P(Y|Z,Theta)) ((9.2))
即：
(P(Y|Theta) = prod_{j=1}^{n}[pi p^{y_j}(1-p)^{1-y_j} + (1-pi)q^{y_j}(1-q)^{1-y_j}]) ((9.3))
考虑求模型参数(Theta=(pi, p, q))的极大似然估计，即：
(hat{Theta} = argmax logP(Y|Theta)) ((9.4))
这个参数只有通过迭代的方法求解。EM算法首先选择参数的初值：
(Theta^{(0)} = (pi^{(0)}, p^{(0)}, q^{(0)}))
然后通过下面的步骤迭代计算参数的估计值，直至收敛。第i次迭代参数的估计值为：
(Theta^{(i)} = (pi^{(i)}, p^{(i)}, q^{(i)}))
则EM算法的第i+1次迭代如下：
E步：计算在模型参数(Theta^{(i)} = (pi^{(i)}, p^{(i)}, q^{(i)}))下观测数据(y_i)来自硬币B的概率：
(mu^{i+1} = pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}/[pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j} + (1-pi^{(i)})(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}]) ((9.5))
M步：计算模型参数的新估计值
(pi^{(i+1)} = 1/nsum_{j=1}^{n}mu_j^{(i+1)}) ((9.6))
(p^{(i+1)} = frac{sum_{j=1}^{n}mu_j^{(i+1)}y_j}{sum_{j=1}^{n}mu_j^{(i+1)}}) ((9.7))
(q^{(i+1)} = frac{sum_{j=1}^{n}(1-mu_j^{(i+1)})y_j}{sum_{j=1}^{n}(1-mu_j^{(i+1)})}) ((9.8))
下面我们一轮一轮来计算：
则每次来自B和C的概率分布如下：
对于硬币B，这10次抛硬币，6次正面，4次反面：p = 0.5 * 6 / (0.5*6 + 0.5 * 4) = 0.6，同理对于硬币C：q = 0.6，而对于A：(pi = 0.5)
下面我们看看第二轮：
则每次来自B和C的概率分布如下：
对于硬币B，这10次抛硬币，6次正面，4次反面：p = 0.5 * 6 / (0.5*6 + 0.5 * 4) = 0.6，同理对于硬币C：q = 0.6，而对于A：(pi = 0.5)
可以看见，已经收敛
下面我们改变初始值：
对于硬币B，这10次抛硬币，6次正面，4次反面：p = 0.3636 * 6 / (0.3636*6 + 0.4706 * 4) = 0.5368，同理对于硬币C：q = 0.6364*6 / (0.6364*6+0.5294*4) = 0.6432
而对于A：(pi =(0.3636*6 + 0.4706*4) / (0.3636*6 + 0.4706*4 + 0.6364*6 + 0.5294*4) = 0.4064 )，即B贡献的正面次数+C贡献的正面次数之和除以10
下面我们看看第二轮：
对于硬币B，这10次抛硬币，6次正面，4次反面：p = 0.3637 * 6 / (0.3637 * 6 + 0.4705 * 4) = 0.5369，同理对于硬币C：q = 0.6363 * 6 / (0.6363*6 + 0.5295*4) = 0.6432
而对于A：(pi =(0.3637*6 + 0.4705*4) / (0.3637*6 + 0.4705*4 + 0.6363*6 + 0.5295*4) = 0.4064 )，即B贡献的正面次数+C贡献的正面次数之和除以10
可以看见，已经收敛
感兴趣的可以参考一下：https://blog.csdn.net/SAJIAHAN/article/details/53106642