EM算法及其推广学习笔记

EM算法及其推广学习笔记

前言

在学习隐马尔科夫模型时，在学习算法中指出了Baum-Welch算法，来实现对隐马尔科夫模型参数的求解。在该学习算法中用到了EM算法，因此我们先来看看EM算法到底是何方神圣。可自己在学习EM算法时，又遇到了一个坑，什么是极大似然函数？因此，本文先介绍极大似然函数的相关概念，然后再对EM算法进行物理映射和实际数学推导。本文需要大量概率论知识，在数学推导关键处会贴出相关参考资料和博文，以备后期查询使用。

思考

1. 什么是极大似然函数，用来解决什么实际问题？
2. EM算法是什么，用来解决什么问题？

正文

极大似然估计

1.概念
在已知试验结果（即是样本）的情况下，用来估计满足这些样本分布的参数，把可能性最大的那个参数$\theta$作为真实$\theta^*$的参数估计。

2.实例
在学习极大似然函数时，一直在思考用什么样的例子才能把概念讲清楚。其实最大的困惑便在于似然估计是假设你有了一堆样本，你需要根据这些样本来猜出某些概率参数，从而使得这些样本在你面前的概率最大。用公式表示为：

$$\max_{\theta}P(X;\theta)$$
这里显然做了一个最基本假设，即出现在我们面前的样本表征了最真实的物理世界，没有之一，且不考虑外界噪声对样本的干扰。现在，咱们举一个实际的例子，来理解似然估计的物理含义。

假设我们已知得肺癌的概率为总体的0.2，而不得肺癌的概率为0.8。在得肺癌中的人群中，我们又做了相应的统计，即抽烟人群占0.7，而不抽烟人群只有0.3。在没有肺癌的人群中，抽烟者占0.2，不抽烟者占0.8。然而很不幸，在真实世界中，我们遇到了一个倒霉蛋，他得了肺癌。问题来了，他是不是抽烟呢？

这个例子不算是真正的似然估计最贴切的实例，但基本思想可以拿来表征一下，方便下面对公式的理解。很显然，得肺癌背后的概率模型显然和抽烟相关性很大，即想要得肺癌的概率最大，模型必须拟合到$\theta=0.7$，使得获肺癌的概率最大。因此，我们可以大胆猜测，该肺癌患者是吸烟人群。详细的关于似然函数的物理含义可以参看博文从最大似然到EM算法浅解。

我们先给出数学定义，极大似然估计可以分为离散型和连续型两类。

3.离散型
设$x$为离散型随机变量，$\theta =(\theta_1,\theta_2,...,\theta_k)$为多维参数向量，如果随机变量$x_1,x_2,...,x_n$相互独立且概率计算式为$P(x=x_i) = P(x_i;\theta_1,...,\theta_k)$，则$P(X=x_1,x_2,...,x_n) = \prod_{i=1}^nP(x_i;\theta_1,...,\theta_k)$.我们就把该式子记作$L(\theta)=\prod_{i=1}^nP(x_i;\theta_1,...,\theta_k)$，称此函数为似然函数。似然函数值的大小意味着该样本值出现的可能性的大小，既然已经得到了样本值$X=(x_1,x_2,...,x_n)$，那么它出现的可能性应该是较大的，即似然函数的值也应该是较大的，因而最大似然估计就是选择使$L(\theta)$达到最大值的那个$\theta$作为真实$\theta^*$的估计。

4.连续型
设$x$为连续型随机变量，概率密度函数为$f(x_i;\theta_1,\theta_2,...,\theta_k),x_1,x_2,...,x_n$为该总体中抽出的样本，同样的如果$x_1,x_2,...,x_n$相互独立且同分布，于是样本的联合概率密度为$L(\theta) = \prod_{i=1}^n f(x_1;\theta_1,\theta_2,...,\theta_k)$，大致过程同离散型一样。

5.应用举例
继续来看例子，假设进行一个实验，实验次数为10次，每次实验成功率为0.2，那么不成功的概率为0.8，用$y$来表示成功的次数。由于前后的实验是相互独立的，所以可以计算得到成功的次数的概率密度为：

$$f(y;\theta = 0.2) = \frac {10!}{y!(10-y)!}(0.2)^y(1-0.2)^{10-y} , y=0,1,...,10.$$

该式子分为两项因子，10次实验中有$y$次成功，那么即在10次中随意挑选$y$个成功的实验，即$C_{10}^y$；第二项为10次实验中，y次实验成功的概率。更一般地，我们可以把每次实验成功的概率当作一个变量$\theta$，则上式可以写为：

$$f(y;\theta) = \frac {10!}{y!(10-y)!}\theta^y(1-\theta)^{10-y} , y=0,1,...,10.$$
显然， $f(y;\theta)$是关于随机变量 $y$和概率参数$\theta$的函数，记作， $L(\theta) = f(y;\theta)$。我们这里由于y是已知变量，所以似然函数只关于参数 $\theta$变化。

当$\theta =0.2$时，我们可以得到y取不同值的概率分布情况，如下图所示：

当$\theta =0.7$时，我们可以得到y取不同值的概率分布情况，如下图所示：

好了，现在假设我们在实验室，开始完成某个实验，我们并不知道该实验成功的概率是多少，但做了10次实验后，我们只成功了2次，用高中的拿点概率知识拿来求解，那不就是实验成功率为0.2。的确，但由于实验次数相当的小，这里的0.2并非是真正的概率，而只是我们实验成功的频率。如抛一枚硬币，抛个10次，可能正面朝上的频率为0.6，但我们都知道，实际正面朝上的概率为0.5。那如何让频率接近0.5呢，不断的增加实验次数即可，你抛个2万次试试。所以我们不能简单的就把这个问题中求解的0.2作为我们的答案，我们也不可能大量重复实验来统计该实验成功率。遇到这种情况，我们便用到了似然估计方法。

如上给出了似然函数：

$$L(\theta)=f(y;\theta) = \frac {10!}{y!(10-y)!}\theta^y(1-\theta)^{10-y} , y=0,1,...,10.$$
现在我们已知实验次数为2，我们要求 $\theta$使得该似然函数取到极大值，求极值问题只需要对 $\theta$求导即可，如果是多参量，那么可以对它求偏导，得到似然方程组，同样能求出想要的解。
$$L'(\theta) = \frac {10!}{2!(10-2)!}\theta^2(1-\theta)^8 = 2-10\theta =0$$
求得 $\theta = 0.2$。咦，算出来的答案是一样的，这不是多此一举嘛，但上述实验成功次数背后的参数 $\theta$模型是一维的，即我们可以用“肉眼”来直接看出答案，假如我们这次不是观察实验成功次数，来猜实验成功率，我换个问题问，假设我们班男生身高符合高斯分布，即男生身高概率密度函数符合高斯分布，给你一群男生的身高，请告诉我高斯分布的方差和均值分别是多少？这种情况下，其背后的 $\theta = [\mu,\pi]$含有两个参数，简单的靠肉眼观察显然无法给出答案，因此我们需要借助数学工具，来理论化的证明说，当看到这一群男生的身高时，我们能找到参数 $\theta= [\mu,\pi]$使得出现 这群男生身高的概率最大，注意是这群而不是那群！言外之意就是说， $\theta$参数能够对男生进行分类！隐约看出了EM算法中的一些思想。

我们再把上述问题复杂一下，假设我们现在重复上述实验过程，即第一次，重复实验10次，观察到实验成功次数为1次；第二次，重复实验10次，观察实验成功次数为2次。问：你能告诉我实验成功的次数为几次吗？还是用数学严格的进行求解一次！

这里我们有两个观察值，即随机变量$y_1=1,y_2=2$，两个随机变量符合相互独立的条件，所以由概率公式得：

$$P(y_1,y_2;\theta) =\frac {10!}{1!(10-1)!}\theta(1-\theta)^9 \times \frac {10!}{2!(10-2)!}\theta^2(1-\theta)^8$$
同样的，要求 $\theta$使得似然函数取极大值，我们需要对等式进行求导，问题来了，这是2个观察值，n个观察值进行求导，那这复杂得根本无法计算。因此简单的想法就是把求导的乘法法则能够映射到求导的加法法则，因此便有了对数似然函数的引出，即取 $log$函数，得：
$$L(\theta) = \log P(X;\theta)$$
这样对上式进行求导便方便很多，更关键的是，求解出来的 $\theta$值与原先的概率分布函数是等价的。
$$\frac {\partial \ln P(y_1,y_2;\theta)}{\partial_\theta} = \frac 1{\theta}-\frac 9{1-\theta} + \frac 2{\theta}-\frac 8{1-\theta} = 0$$
求得 $\theta = 0.15$。即试验成功的概率为0.15。

6.总结
最大似然估计，只是一种概率论在统计学的应用，它是参数估计的方法之一。说的是已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，参数估计就是通过若干次实验，观察其结果，利用结果推出参数的大概值。对似然的理解是基于最基本的假设，我们得到的观察结果是在某个参数模型下出现概率最大的。我们并不考虑，实际当实验成功率为0.7时，我们观察两次，分别只观察到10次实验中成功1次和2次。这种情况在某种环境下可能会发生，但倘若发生，我们就认为参数$\theta =0.15$，而不是0.7。
求最大似然函数估计值的一般步骤：
1. 写出似然函数
2. 对似然函数取对数，并整理
3. 求导数
4. 解似然方程

上述内容摘自博文–最大似然估计学习总结——MadTurtle，有省略的部分也有补充的部分。

EM算法

1.定义
概率模型有时即含有观测变量，又含有隐变量或潜在变量。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。我们仅讨论极大似然估计，极大后验概率估计与其类似。（隐含了一个问题，概率模型中当存在隐变量时，就无法直接用极大似然估计法进行求解，这是为什么？）

2.三硬币模型
假设有3枚硬币，分别记作A，B，C。这些硬币正面出现的概率分别是$\pi,p,q$.进行如下投掷试验：先掷硬币A，根据其结果选出硬币B或硬币C，正面选硬币B，反而选硬币C；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立地重复n次试验（这里，n=10），预测结果如下：

$$1,1,0,1,0,0,1,0,1,1$$
假设智能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数。

同样的，先用先前似然估计方法来求解一波，看看能否给出答案。假设我们知道了一个观测值：

$$P(y|\theta) = \sum_z P(y,z|\theta) = \sum_z P(z|\theta)P(y|z,\theta) \space 取y =0,1$$
直接用 $\theta = (\pi,p,q)$得
$$P(y|\theta) = \pi p^y(1-p)^{1-y} + (1-\pi) q^y(1-q)^{1-y}$$
该式子中 $y$为已知，其他参数均位置，假设我们知道观察序列的第一次投掷结果为1，因此把$y =1$代入得
$$P(y|\theta) = \pi p + (1-\pi)q$$
以极大似然方法进行求解，分别对参数 $\pi,p,q$进行求导，你会发现对 $\pi$求导，求出 $p=q$来，对 $p,q$分别求导求出 $\pi =0 和 \pi =1$，显然是没有解析解。因此，传统的似然估计方法是无法解决上述三硬币模型的问题。这也是为什么EM算法提出的原因，即它能解决传统求导解决不了的问题。（遗留一个问题和一个证明，三硬币模型中是由于 $\pi$的出现，使得似然方程无解？是一旦隐藏变量出现，就无法求解似然方程嘛？如何证明。）

这里，随机变量$y$是观测变量，表示一次试验观测的结果是1或0；随机变量$z$是隐变量，表示未观测到的掷硬币A的结果；$\theta =(\pi,p,q)$是模型参数。这一模型是以上数据的生成模型。注意，随机变量$y$的数据可以预测，随机变量$z$的数据不可观测。（我们无法得知掷硬币A是正面还是反面，信息缺失，数学建模，无法求出解析解？）

将观测数据表示为 $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)$$
即
$$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}]$$
考虑求模型参数 $\theta = (\pi,p,q)$的极大似然估计，即
$$\hat \theta = arg \max_{\theta} \log P(Y|\theta)$$
$\log$函数中有加法，对参数进行求偏导显然是困难的。因此，我们需要另辟蹊径来求解该似然函数的极大值。也就是该似然函数并非是单纯的观测随机变量的概率分布函数，而是隐藏了不可观测变量的概率分布，显然如果并不知道每个样本背后隐藏变量的值，那么求解出来的参数是无意义的。

回到身高分布问题，我们现在假设在全校抽中了100个男生，100个女生，男生和女生的身高都符合正态分布，即背后的概率模型为正态分布概率密度函数，对样本进行统计时，我们标注了他是男生，她是女生，且观测到了每个人的身高，据此我们就可以列出一个关于男生和女生的似然函数，根据似然估计方法，我们通过求导的方式便能求解出模型的参数$\theta = (\mu,\delta)^T$。

但假设我们在全校抽中了200个人，但背后我们却没有统计它们的性别，如果把这200个人当作同一群人进行似然函数的建模求解，那显然是不太明智的，因为我们都知道，性别的差异对身高的影响是相当大的，我们非要把女生当男生，并由此计算出模型的参数来，再通过该模型预测出来的男生身高，就显得不那么准确了。所以，在针对隐藏变量的情况下，我们需要考虑性别这样的因素，这也就有了EM算法提出的意义。

在对问题进行分析时，其实我们是可以不用考虑隐藏变量的因素的？但诸如隐马尔可夫模型中的学习问题时，很显然它有隐含的状态，那么进行参数学习时，我们就需要用到EM算法。同样地，三硬币模型也明确的告诉我们含有隐藏变量，因此也必须使用EM算法进行求解。

刚才说了，对含有隐藏变量的似然函数是无法用求导的方式进行求解的，我们先把式子写出来，即

$$L(\theta) = \sum_i \log p(x^{(i)};\theta) = \sum_i \log \sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)$$
EM算法是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化 $L(\theta)$，我们可以不断地建立 $l$的下界（E步），然后优化下界（M步）。

对于每一个样例$i$，让 $Q_i$表示该样例隐含变量 $z$的某种分布，$Q_i$满足的条件是 $\sum_zQ_i(z)=1,Q_i(z) \ge 0.$比如要将学校抽的200人进行聚类，假设隐藏变量 $z$是身高，那么就是连续的高斯分布。如果按照隐藏变量是男女，那么就是伯努利分布。可以有前面阐述的内容得到下面的公式：

$$\sum_i \log p(x^{(i)};\theta) = \sum_i \log \sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta) \\ =\sum_i \log \sum_{z^{(i)}}Q_i(z^{(i)}) \frac {p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \\ \ge \sum_i \sum_{z^{(i)}}Q_i(z^{(i)})\log \frac {p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$$
第一步和第二步比较直接，就是分子分母同乘以一个相等的函数。第二步和第三步利用了Jesson不等式，具体的推导过程请参看博文- EM算法原理。

这个过程可以看作是对$L(\theta)$求了下界。对于$Q_i$的选择，有多种可能，哪种更好呢？假设$\theta$给定，那么$L(\theta)$的值就决定于$Q_i(z^{(i)}$和$p(x^{(i)},z^{(i)};\theta)$。我们可以通过调整这两个概率是下界不断上升，以逼近$L(\theta)$的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率都能够等价于$L(\theta)$了。按照这个思路，我们要找到等式成立的条件。根据Jessen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：

$$\frac {p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} =C$$
C为常数，不依赖于 $z^{(i)}$。对此式子做进一步推导，我们知道 $\sum_ZQ_i(z^{(i)}) =1$，可以推导得：
$$Q_i(z^{(i)}) = p(z^{(i)}|x^{(i)}; \theta)$$
对推导过程感兴趣的，可以继续参看博文- EM算法原理。

这里简单提一下取常数的物理含义，上图即为 $\log$函数，而在 $\log$函数上两个点，可以分别用 $\frac {p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$和 $\frac {p(x^{(j)},z^{(j)};\theta)}{Q_j(z^{(j)})},i,j$分别表示不同的样本，要让Jessen不等式成立，显然需要让这两点重合，才能取得等号。具体的解释请参看视频教程 七月在线-18分钟理解EM算法。

刚开始接触这个$L(\theta)$不理解为什么要让下界函数和$L(\theta)$等号成立，然后再开始对新的下界函数进行求极值的过程，并且求极值过程能够逼近$L(\theta)$的极值，不急，咱们来看看博文从最大似然到EM算法浅解的解释，参看下图：

E步的过程，就是调整$Q(z)$使得下界$J(Z,Q)$不断上升，直到与$L(\theta)$在$\theta$点重合，找到了下界函数后，固定$Q(z)$，对参数进行迭代，找到$J(Z,\theta)$的极大值，反复上述操作，从而逼近$L(\theta)$的极大值。（两个问题，$\theta$怎么变？该过程为何收敛？）参看书本《统计学习方法》第159页

算法（EM算法）

输入：观测变量数据Y，隐变量数据Z，联合分布$P(Y,Z|\theta)$，条件分布$P(Z|Y,\theta)$
输出：模型参数$\theta$
(1) 选择参数的初值$\theta^{(0)}$，开始迭代；
(2) E步：记$\theta^{(0)}$为第$i$次迭代参数$\theta$的估计值，在第$i+1$次迭代的E步，计算

$$Q(\theta,\theta^{(i)}) = \sum_Z \log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})$$
这里， $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_{\theta} Q(\theta,\theta^{(i)})$$
(4) 重复第(2)步和第(3)步，直到收敛。

Q函数的推导在书本《统计学习方法》第159页。

行文至此，除了在数学上能够强行“理解”证明的过程，但实际的物理问题却始终无法对应到这些形式化的符号中去，暂且把这些恼人数学放一边，来看看针对某些特定的实际问题，算法是如何解决的，没准能够从中理解一些数学公式的实际含义。

Code Time

双硬币模型

假设有两枚硬币A、B，以相同的概率随机选择一个硬币，进行如下的抛硬币实验：共做5次实验，每次实验独立的抛十次，结果如图中a所示，例如某次实验产生了H、T、T、T、H、H、T、H、T、H，H代表正面朝上。

假设试验数据记录员可能是实习生，业务不一定熟悉，造成a和b两种情况

a表示实习生记录了详细的试验数据，我们可以观测到试验数据中每次选择的是A还是B

b表示实习生忘了记录每次试验选择的是A还是B，我们无法观测实验数据中选择的硬币是哪个

问在两种情况下分别如何估计两个硬币正面出现的概率？

这是实习生a记录的情况，由于这里数据的背后参考模型已知（已分好类），因此用极大似然估计方法就能分别求出$\hat \theta_A和\hat \theta_B$的概率来。与上文第一节中的例子完全类似。

这是实习生b记录的情况，令人糟糕的是数据的背后参考模型混淆在了一起，我们无法得知它们这几组实验数据是由A抛的还是由B抛的，因为这里隐含了一个该组数据是A还是B的分类问题。抛开数学对隐含变量的复杂求解过程，我们可以先给出一个思路来解决上述问题。

第一，既然我们不清楚是A类还是B类，但假定我们初始化了A类硬币抛正面的概率和B类硬币抛正面的概率，这两者的概率是随意定的，但由于两者概率可以存在差异，假设$P(y = H;\theta_A) \gt P(y = H;\theta_B)$，那么一个明显的特征就是，由于能观察到10次硬币中有几次是成功的，我们可以基于这次观察，给出$P(z = A| Y;\theta_A,\theta_B)$的概率，上式的含义是可以根据两个参数的初值求出，在给定观察序列的情况下，它属于A类还是B类的概率。用公式可以表示为：

$$P(z | Y;\theta_A,\theta_B) = \frac {P(Y,z;\theta_A,\theta_B)}{P(Y;\theta_A,\theta_B)}$$
其中，z表示单个观察到的随机变量，此处z=A or B（属于分类问题），Y表示观察序列，即 $Y =(y_1,y_2,...,y_{10})^T$。由此，给定观察序列后，我们可以算出属于A类的概率和属于B类的概率，那么很显然CoinA 和CoinB 不再是属于你有我没有，你没有我有的敌对关系，因为我自己都还不是很清楚是不是A类，由此10个硬币，我们就根据概率进行一次平均分配咯，这样CoinA 和CoinB 在一次观察结果时，都能得到属于自己的那一份，非常的和谐。这一部分便是求期望的过程，即对于第一个观察序列中，10次抛硬币过程中5次为正面朝上，令 $y_j = 5$，由此可以得到关于隐含变量的数学期望 $E(z) = 0.45* 5 +0.55 * 5$，“+”号的左边右边，分别表示CoinA的分配和CoinB的分配。分配的份额根据z函数的分布给定，z函数的分布规则根据缺失的信息进行建模，解由初始参数求出。

因此分类问题，给出了每个CoinA 和CoinB 的分配额，有了所有观察值CoinA和CoinB的分配额，我们就可以单独的对CoinA和CoinB进行最大似然估计方法。求出来的新参数，再求z函数，求期望，求参数，如此迭代下去，直到收敛到某一结果。

算法实现

在双硬币模型中，对某个种类硬币投掷10次中成功n次概率模型$P(y|\theta)= \binom{10}{n}\theta^n(1-\theta)^{(10-n)},y =0,1,...,10$符合伯努利分布。

前期准备工作
我们可以实际的操作一把，来看看在python中如何可视化伯努利分布。前期环境搭建：
1. 安装python扩展包：scipy,matplotlib,seaborn,ipywidgets。
- 在安装扩展包过程中遇到了一些课，如安装scipy，用pip直接安装遇到安装失败的情况，如果有类似的情况可以参考链接[Python]Windows7 x64安装numpy和scipy
2. 安装ipython，貌似在python3.5版本中，自带了ipython。
- 安装完毕，在cmd中敲入ipython notebook，开始编程之旅

进入游览器

点击new按钮后，新增python3按钮，进入下一界面

实际编写
接下类就可以进行我们实际的操作了，首先引入一些必须的科学计算包：

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
from ipywidgets import interact, FloatSlider
%matplotlib inline

编写伯努利分布函数

a=range(11)
def plot_binomial(p=0.5):
    fig, ax = plt.subplots(figsize=(4,3))
    y = [0]*11
    for i in a:
        y[i-1] =  stats.binom.pmf(i, 10, p)
    ax.bar(a,y,label="$p = %.1f$" % p)
    ax.set_ylabel('PMF of $k$ heads')
    ax.set_xlabel('$k$')
    ax.set_ylim((0,0.5))
    ax.set_xlim((0,10))
    ax.legend()
    return fig

可视化伯努利分布函数

interact(plot_binomial, p=FloatSlider(min=0.0,max=1.0,step=0.1,value=0.2))

可视化界面如下

这是在p=0.2的情况下的伯努利分布函数，代回双硬币模型中去，当观察到10次实验中只有2次成功了，那么该$\theta$参数便是0.2。因为只有当$\theta =0.2$时，10次实验中出现成功次数为2次的概率最大。

回到实习生b，由数据可得观察矩阵为：

observations = np.array([[1,0,0,0,1,1,0,1,0,1],
                         [1,1,1,1,0,1,1,1,1,1],
                         [1,0,1,1,1,1,1,0,1,1],
                         [1,0,1,0,0,0,1,1,0,0],
                         [0,1,1,1,0,1,1,1,0,1]])

实际每组观察数据属于A类，B类的隐藏状态为：

coins_id = np.array([False,True,True,False,True])

那么在观察数组中，属于A类的数组为：

observations[coins_id]

输出结果为：

array([[1, 1, 1, 1, 0, 1, 1, 1, 1, 1],
       [1, 0, 1, 1, 1, 1, 1, 0, 1, 1],
       [0, 1, 1, 1, 0, 1, 1, 1, 0, 1]])

在所有属于A类的数组中，总的实验成功次数为：

np.sum(observations[coins_id])

输出结果为：

24

所以说，针对属于A类的硬币，它的参数$\theta_A$：

1.0*np.sum(observations[coins_id])/observations[coins_id].size

输出结果为：

0.80000000000000004

同理，对于属于B类的硬币，它的参数为$\theta_B$：

1.0*np.sum(observations[~coins_id])/observations[~coins_id].size

输出结果为：

0.45000000000000001

EM算法步骤

首先来看看，针对第一组观察数据，每一步的数据是如何求出的。

# 对于实验成功率为0.6的情况，10次中成功5次的概率
coin_A_pmf_observation_1= stats.binom.pmf(5,10,0.6)
coin_A_pmf_observation_1

# 均是针对第一组观察数据的情况
coin_B_pmf_observation_1= stats.binom.pmf(5,10,0.5)
coin_B_pmf_observation_1

# 针对第一组观察数据中，属于A类硬币投掷的概率p(z=A|Y,theta)
normalized_coin_A_pmf_observation_1 = coin_A_pmf_observation_1/(coin_A_pmf_observation_1+coin_B_pmf_observation_1)
print ("%0.2f" %normalized_coin_A_pmf_observation_1)

# 针对第一组观察数据中，属于B类硬币投掷的概率p(z=B|Y,theta)
normalized_coin_B_pmf_observation_1 = coin_B_pmf_observation_1/(coin_A_pmf_observation_1+coin_B_pmf_observation_1)
print ("%0.2f" %normalized_coin_B_pmf_observation_1)

#求期望过程
weighted_heads_A_obervation_1 = 5*normalized_coin_A_pmf_observation_1
print ("Coin A Weighted count for heads in observation 1: %0.2f" %weighted_heads_A_obervation_1)
weighted_tails_A_obervation_1 = 5*(1-normalized_coin_A_pmf_observation_1)
print ("Coin A Weighted count for tails in observation 1: %0.2f" %weighted_tails_A_obervation_1)
weighted_heads_B_obervation_1 = 5*normalized_coin_B_pmf_observation_1
print ("Coin B Weighted count for heads in observation 1: %0.2f" %weighted_heads_B_obervation_1)
weighted_tails_B_obervation_1 = 5*(1-normalized_coin_B_pmf_observation_1)
print ("Coin B Weighted count for tails in observation 1: %0.2f" %weighted_tails_B_obervation_1)

单步EM算法，迭代一次算法实现步骤

def em_single(priors,observations):
    """
    performs a single Em step
    Arguments
    ---------
    priors:[theta_A,theta_B]
    observations:[m X n matrix]

    Returns
    -------
    new_priors:[new_theta_A,new_theta_B]
    """
    counts ={'A':{'H':0,'T':0},'B':{'H':0,'T':0}}
    theta_A = priors[0]
    theta_B = priors[1]

    # E setp
    for observation in observations:
        len_observation = len(observation)
        # 对数组求和 head =1，tail =0
        num_heads = observation.sum()
        num_tails = len_observation - num_heads
        # 如果属于A类，A的分布概率 和属于B类的分布概率
        contribution_A = stats.binom.pmf(num_heads,len_observation,theta_A)
        contribution_B = stats.binom.pmf(num_heads,len_observation,theta_B)
        # z分布概率
        weight_A = contribution_A /(contribution_A+contribution_B)
        weight_B = contribution_B /(contribution_A+contribution_B)

        # Incrementing counts
        counts['A']['H'] += weight_A*num_heads
        counts['A']['T'] += weight_A*num_tails
        counts['B']['H'] += weight_B*num_heads
        counts['B']['T'] += weight_B*num_tails

    # M setp
    new_theta_A = counts['A']['H']/(counts['A']['H']+counts['A']['T'])
    new_theta_B = counts['B']['H']/(counts['B']['H']+counts['B']['T'])

    return [new_theta_A, new_theta_B]

迭代一次输出结果为：

em_single([0.6,0.5],observations)
# 结果
[0.71301223540051617, 0.58133930831366265]

收敛算法

def em(observations,prior,tol=1e-6,iterations=1000):
    import math
    iteration = 0
    while iteration <iterations:
        new_prior = em_single(prior,observations)
        delta_change = np.abs(prior[0] -new_prior[0])
        if delta_change < tol:
            break
        else:
            prior = new_prior
            iteration +=1
    return [new_prior,iteration]

最终结果为：

em(observations, [0.6,0.5])
# 结果
[[0.79678875938310978, 0.51958393567528027], 14]

详细代码请参看博文Programatically understanding Expectation Maximization algorithm

最终实习生b的EM算法得出的结果，跟实习生a得出的参数还是非常相近的，但EM算法跟初始值的设置有着很大的关系，不信，修改[06,0.5]多尝试尝试。

针对EM算法的所有相关内容都已经阐述完毕了，在学习过程中，遇到了许多坑，诸如数学公式与实际的问题没法完全一一对应，但个人认为，自家数学功底不够深厚时，唯有参照实际情况，来解决一个实际问题，多加练习，在日后回过头来看数学定义时，没准能够恍然大悟。解决双硬币模型的思想是很简单的，但对当细节问题进行深入理解时，本文还有很多不足，如算法如何收敛，为何收敛，算法E-step为何是求期望。待补足数学上的漏洞或者经历过大量实践后，有望能够解决，未完待续。

参考文献

• 从最大似然到EM算法浅解
• 最大似然估计学习总结——MadTurtle
• EM算法原理
• 七月在线-18分钟理解EM算法
• EM算法及其推广 - 码农场