<强>前言:强>前一篇文章大概说了EM算法的整个理解以及一些相关的公式神马的,那些数学公式啥的看完真的是忘完了,那就来用代码记忆记忆吧!接下来将会对python版本的EM算法进行一些分析。
<强> EM的python实现和解析强>
<强>引入问题(双硬币问题)强>
假设有两枚硬币A, B,以相同的概率随机选择一个硬币,进行如下的抛硬币实验:共做5次实验,每次实验独立的抛十次,结果如图中一个所示,例如某次实验产生了H, T, T, T, H, H, T, H, T, H, H代表正面朝上。
假设试验数据记录员可能是实习生,业务不一定熟悉,造成a和b两种情况
表示实习生记录了详细的试验数据,我们可以观测到试验数据中每次选择的是一个还是B
b表示实习生忘了记录每次试验选择的是一个还是b,我们无法观测实验数据中选择的硬币是哪个
问在两种情况下分别如何估计两个硬币正面出现的概率?
以上的针对于b实习生的问题其实和三硬币问题类似,只是这里把三硬币中第一个抛硬币的选择换成了实习生的选择。
对于已知是一个硬币还是B硬币抛出的结果的时候,可以直接采用概率的求法来进行求解。对于含有隐变量的情况,也就是不知道到底是一个硬币抛出的结果还是B硬币抛出的结果的时候,就需要采用EM算法进行求解了。如下图:
其中的EM算法的第一步就是初始化的过程,然后根据这个参数得出应该产生的结果。
<>强构建观测数据集强>
针对这个问题,首先采集数据,用1表示H(正面),0表示T(反面):
#硬币投掷结果 观察=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行数据的μ。
有了μ,就可以估计数据中AB分别产生正反面的次数了.μ代表数据来自硬币一个的概率的估计,将它乘上正面的总数,得到正面来自硬币一个的总数,同理有反面,同理有B的正反面。
#更新在当前参数下,B硬币产生的正反面次数 数[A] [' H '] +=weight_A * num_heads 数[A] [' T '] +=weight_A * num_tails 数[B] [' H '] +=weight_B * num_heads 数[B] [' T '] +=weight_B * num_tails<强>第一个迭代的M步强>
当前模型参数下,AB分别产生正反面的次数估计出来了,就可以计算新的模型参数了:
new_theta_A=计数[A] [' H ']/(计数[A] [' H '] +计数[A] [' T ']) new_theta_B=计数[B] [' H ']/(计数[B] [' H '] +计数[B] [' T '])于是就可以整理一下,给出EM算法单个迭代的代码:
def em_single(先知先觉,观察): ”“” 他们算法的单次迭代 参数 ------------ 先验(theta_A, theta_B): 观察(m X n矩阵): 返回 --------------- new_priors (new_theta_A, new_theta_B): :param先验: :param观察: 返回: ”“” 数量={A: {“H”: 0,“T”: 0}, B: {“H”: 0,“T”: 0}} theta_A=先验[0] theta_B=先验[1] # E步骤 在观察观察: len_observation=len(观察) num_heads=observation.sum () num_tails=len_observation-num_heads #二项分布求解公式 contribution_A=scipy.stats.binom.pmf (num_heads len_observation theta_A) contribution_B=scipy.stats.binom.pmf (num_heads len_observation theta_B) weight_A=contribution_A/(contribution_A + contribution_B) weight_B=contribution_B/(contribution_A + contribution_B) #更新在当前参数下,B硬币产生的正反面次数 数[A] [' H '] +=weight_A * num_heads 数[A] [' T '] +=weight_A * num_tails 数[B] [' H '] +=weight_B * num_heads 数[B] [' T '] +=weight_B * num_tails #米一步 new_theta_A=计数[A] [' H ']/(计数[A] [' H '] +计数[A] [' T ']) new_theta_B=计数[B] [' H ']/(计数[B] [' H '] +计数[B] [' T ']) 返回(new_theta_A, new_theta_B)他们算法的python实现的方法步骤