EM算法的python实现的方法步骤

前言：前一篇文章大概说了EM算法的整个理解以及一些相关的公式神马的，那些数学公式啥的看完真的是忘完了，那就来用代码记忆记忆吧！接下来将会对python版本的EM算法进行一些分析。

EM的python实现和解析

引入问题（双硬币问题）

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

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

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

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

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

以上的针对于b实习生的问题其实和三硬币问题类似，只是这里把三硬币中第一个抛硬币的选择换成了实习生的选择。

对于已知是A硬币还是B硬币抛出的结果的时候，可以直接采用概率的求法来进行求解。对于含有隐变量的情况，也就是不知道到底是A硬币抛出的结果还是B硬币抛出的结果的时候，就需要采用EM算法进行求解了。如下图：

其中的EM算法的第一步就是初始化的过程，然后根据这个参数得出应该产生的结果。

构建观测数据集

针对这个问题，首先采集数据，用1表示H（正面），0表示T（反面）：

#硬币投掷结果observations = numpy.array([[1,0,0,0,1,1,0,1,0,1],            [1,1,1,1,0,1,1,1,0,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]])

第一步：参数的初始化

参数赋初值

第一个迭代的E步

抛硬币是一个二项分布，可以用scipy中的binom来计算。对于第一行数据，正反面各有5次，所以：

#二项分布求解公式contribution_A = scipy.stats.binom.pmf(num_heads,len_observation,theta_A)contribution_B = scipy.stats.binom.pmf(num_heads,len_observation,theta_B)

将两个概率正规化，得到数据来自硬币A，B的概率：

weight_A = contribution_A / (contribution_A + contribution_B)weight_B = contribution_B / (contribution_A + contribution_B)

这个值类似于三硬币模型中的μ，只不过多了一个下标，代表是第几行数据（数据集由5行构成）。同理，可以算出剩下的4行数据的μ。