EM学习-思想和代码

EM算法的简明实现

当然是教学用的简明实现了，这份实现是针对双硬币模型的。

双硬币模型

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

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

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

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

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

a情况相信大家都很熟悉，既然能观测到试验数据是哪枚硬币产生的，就可以统计正反面的出现次数，直接利用最大似然估计即可。

b情况就无法直接进行最大似然估计了，只能用EM算法，接下来引用nipunbatra博主的简明EM算法Python实现。

 # -*- coding: utf-8 -*-

 """

 Created on Tue Jul  4 18:23:28 2017

 @author: Administrator

 """

 import numpy as np

 from scipy import stats

 priors = [0.6, 0.5]

 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]])

 def em_single(priors, observations):

     """

     input:

         priors:[theta_A, theta_B]

         obvervations:m*n matrix

     output:

     """

     theta_A = priors[0]

     theta_B = priors[1]

     counts = {'A':{'H':0,'T':0}, 'B':{'H':0,'T':0}}

     # e-step

     for observation in observations:

         len_observation = len(observation)

         num_heads = observation.sum()   # 正面个数

         num_tails = len_observation - num_heads     # 反面个数

         # 两个二项分布

         contribution_A = stats.binom.pmf(num_heads, len_observation, theta_A)

         contribution_B = 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)

         # 更新在当前参数下，硬币A和B产生正反面的次数

         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-step

     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]

 def em(observations, prior, tol=1e-6, iterations=10000):

     """

     EM算法

     param observations: 观察数据

     param prior: 模型初值

     param tol: 迭代结束阈值

     param iteration: 最大迭代数

     return: 局部最优的模型参数

     """

     import math

     iter = 0

     while iter < iterations:

         new_prior = em_single(prior, observations)

         delta_change = np.abs(new_prior[0]-prior[0])

         if delta_change < tol:

             break

         else:

             prior = new_prior

             iter += 1

         print (iter)

     return [new_prior, iter]

 y = em(observations, priors)

参考自：http://www.hankcs.com/ml/em-algorithm-and-its-generalization.html

巴特西

EM学习-思想和代码

EM算法的简明实现

双硬币模型

最新文章

热门文章