EM算法双硬币模型的python实现
1 双硬币模型
$\quad`假设有两枚硬币A、B,以相同的概率随机选择一个硬币,进行如下的抛硬币实验:共做5次实验,每次实验独立的抛10次,结果如图中a所示,例如某次实验产生了H、T、T、T、H、H、T、H、T、H,H代表正面朝上。
假设试验数据记录员可能是实习生,业务不一定熟悉,造成下面两种情况 :
$\quad`a) 表示实习生记录了详细的试验数据,我们可以观测到试验数据中每次选择的是A还是B 。
$\quad`b) 表示实习生忘了记录每次试验选择的是A还是B,我们无法观测实验数据中选择的硬币是哪个 。
KaTeX parse error: Can't use function '$' in math mode at position 6: \quad$̲`问在两种情况下分别如何估计两…\quad$`以上的针对于b)实习生的问题其实和三硬币问题类似,只是这里把三硬币中第一个抛硬币的选择换成了实习生的选择。
$\quad`对于已知是A硬币还是B硬币抛出的结果的时候,可以直接采用概率的求法来进行求解。对于含有隐变量的情况,也就是不知道到底是A硬币抛出的结果还是B硬币抛出的结果的时候,就需要采用EM算法进行求解了。如下图:
2 python实现
(1) 构建观测数据集
$\quad`针对这个问题,首先采集数据,用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]])
(2) 第一步:参数的初始化
theta_A=0.6;theta_B=0.5
(3) 第一个迭代的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)
有了weight_A,weight_B,就可以估计数据中A、B分别产生正反面的次数了。 weight_A 代表数据来自硬币A的概率的估计,将它乘上正面的总数,得到正面来自硬币A的总数,同理有反面,同理有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
(4) 第一个迭代的M步
当前模型参数下,A、B分别产生正反面的次数估计出来了,就可以计算新的模型参数了:
new_theta_A = counts['A']['H']/(counts['A']['H'] + counts['A']['T'])
new_theta_B = counts['B']['H']/(counts['B']['H'] + counts['B']['T'])
于是就可以整理一下,给出EM算法单个迭代的代码:
def em_single(priors,observations):""" EM算法的单次迭代 Arguments ------------ priors:[theta_A,theta_B] observation:[m X n matrix] Returns --------------- new_priors:[new_theta_A,new_theta_B] :param priors: :param observations: :return: """counts = {'A': {'H': 0, 'T': 0}, 'B': {'H': 0, 'T': 0}}theta_A = priors[0]theta_B = priors[1]# E stepfor observation in observations:len_observation = len(observation)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)# 更新在当前参数下A,B硬币产生的正反面次数counts['A']['H'] += weight_A * num_headscounts['A']['T'] += weight_A * num_tailscounts['B']['H'] += weight_B * num_headscounts['B']['T'] += weight_B * num_tails# M stepnew_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]
(5) EM算法主循环
给定循环的两个终止条件:模型参数变化小于阈值;循环达到最大次数,就可以写出EM算法的主循环了.
def em(observations,prior,tol = 1e-6,iterations=10000):""" EM算法 :param observations:观测数据 :param prior:模型初值 :param tol:迭代结束阈值 :param iterations:最大迭代次数 :return:局部最优的模型参数 """iteration = 0;while iteration < iterations:new_prior = em_single(prior,observations)delta_change = numpy.abs(prior[0] - new_prior[0])if delta_change < tol:breakelse:prior = new_prioriteration += 1return [new_prior,iteration]
(6) 调用及结果
给定数据集和初值,就可以调用EM算法了:
import numpy
import scipy
import scipy.stats
# 硬币投掷结果
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]])
em(observations,[0.6,0.5])
得到:
[[0.72225028549925996, 0.55543808993848298], 36]
我们可以改变初值,试验初值对EM算法的影响。
em(observations,[0.5,0.6])
结果为:
[[0.55543727869042425, 0.72225099139214621], 37]
看来EM算法还是很健壮的。如果把初值设为相等会怎样?
em(observations,[0.3,0.3])
输出:
[[0.64000000000000001, 0.64000000000000001], 1]
显然,两个值相加不为1的时候就会破坏这个EM函数。
换一下初值:
em(observations,[0.99999,0.00001])
输出:
[[0.72225606292866507, 0.55543145006184214], 33]
EM算法对于参数的改变还是有一定的健壮性的。
EM算法双硬币模型的python实现相关推荐
- 机器学习算法_机器学习之EM算法和概率图模型
[晓白]今天我准备更新Machine Learning系列文章希望对机器学习复习和准备面试的同学有帮助!之前更新了感知机和SVM,决策树&代码实战,关注我的专栏可以的文章哦!今天继续更新EM算 ...
- EM 算法与 GMM 模型
EM算法与GMM模型 – 潘登同学的Machine Learning笔记 文章目录 EM算法与GMM模型 -- 潘登同学的Machine Learning笔记 GMM模型 单高斯模型 GM的参数估计( ...
- 机器学习之EM算法的原理及推导(三硬币模型)及Python实现
EM算法的简介 EM算法由两步组成:E步和M步,是最常用的迭代算法. 本文主要参考了李航博士的<统计学习方法> 在此基础上主要依据EM算法原理补充了三硬币模型的推导. 1.EM算法的原理 ...
- 11.EM算法、HMM模型
EM算法入门 算法介绍 极⼤似然估计 EM算法实例描述 EM算法流程 EM算法实例 EM初级版 EM进阶版 HMM模型入门 马尔科夫链 ⻢尔科夫链即为状态空间中从⼀个状态到另⼀个状态转换的随机过程. ...
- baum welch java_Baum-Welch算法(EM算法)对HMM模型的训练
Baum-Welch算法就是EM算法,所以首先给出EM算法的Q函数 \[\sum_zP(Z|Y,\theta')\log P(Y,Z|\theta) \] 换成HMM里面的记号便于理解 \[Q(\la ...
- 数学建模算法(1)—规划模型及其python实现
规划问题及其Python求解方法 什么是规划问题 规划问题的分类 线性规划及其Python解法 非线性规划及其Python解法 什么是规划问题 在人们的生产实践中,经常会遇到如何利用现有资源来安排生产 ...
- EM算法初探——公式推导和三硬币模型解析
EM算法初探--公式推导和三硬币模型解析 转载借鉴:http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html#!comments ...
- 机器学习最易懂之EM算法详解与python实现
文章目录 0.前言 1.EM算法引入 2.具体的EM算法 3.EM算法推导 3.1 Jensen不等式 3.2 EM推导 3.3 EM算法的收敛性 4.EM算法在高斯混合模型中的应用 4.1 高斯混合 ...
- EM算法-硬币实验的理解
EM算法-使用硬币实验的例子理解 EM算法,即最大期望算法(Expectation-Maximization algorithm, EM),是一类通过迭代进行极大似然估计(Maximum Likeli ...
最新文章
- 37 函数的定义和调用
- OSM OpenStreetMap 获取城市路网数据及转为ESRI shp数据的方法
- Matlab读取点云数据显示
- HDU Senior's Gun (水题)
- SSM项目使用GoEasy 实现web消息推送服务
- flink的datastream进行join操作没有输出结果一例
- C语言变量定义和赋值
- 51位院士同写一本书——《两院院士忆高考》新书发布
- 原码、反码、补码,计算机中负数的表示
- 研究员发现70个web缓存投毒漏洞,获奖4万美元
- ocm认证年薪多少_从复读才考上三本,到华为201万年薪的天才少年,他经历了什么?...
- [读书笔记]机器学习:实用案例解析(4)
- python学习第21天
- vector,list,deque容器的迭代器简单介绍
- atitit.web 推送实现方案集合(2)---百度云,jpush 极光推送 ,个推的选型比较.o99
- Bandicam v5.2.1.1860 班迪录屏绿色便携版
- android 5播放flash插件下载地址,Flash Player安卓版
- Qt:Exception at 0xeefde9, code:0x0000005: read access violation at: 0x0, flags = 0x0(first chance)
- AI面试锦囊|网易互娱AI Lab人工智能研究工程师两面分享
- HTML5期末大作业:旅游网页设计与实现——四川成都-(9页 带购物车)