Python实现EM算法实例代码
作者:程序员大本营 发布时间:2021-05-06 03:02:26
EM算法实例
通过实例可以快速了解EM算法的基本思想,具体推导请点文末链接。图a是让我们预热的,图b是EM算法的实例。
这是一个抛硬币的例子,H表示正面向上,T表示反面向上,参数θ表示正面朝上的概率。硬币有两个,A和B,硬币是有偏的。本次实验总共做了5组,每组随机选一个硬币,连续抛10次。如果知道每次抛的是哪个硬币,那么计算参数θ就非常简单了,如
下图所示:
如果不知道每次抛的是哪个硬币呢?那么,我们就需要用EM算法,基本步骤为:
1、给θ_AθA和θ_BθB一个初始值;
2、(E-step)估计每组实验是硬币A的概率(本组实验是硬币B的概率=1-本组实验是硬币A的概率)。分别计算每组实验中,选择A硬币且正面朝上次数的期望值,选择B硬币且正面朝上次数的期望值;
3、(M-step)利用第三步求得的期望值重新计算θ_AθA和θ_BθB;
4、当迭代到一定次数,或者算法收敛到一定精度,结束算法,否则,回到第2步。
计算过程详解:初始值θ_A^{(0)}θA(0)=0.6,θ_B^{(0)}θB(0)=0.5。
由两个硬币的初始值0.6和0.5,容易得出投掷出5正5反的概率是p_A=C^5_{10}*(0.6^5)*(0.4^5)pA=C105∗(0.65)∗(0.45),p_B=C_{10}^5*(0.5^5)*(0.5^5)pB=C105∗(0.55)∗(0.55), p_ApA/(p_ApA+p_BpB)=0.449, 0.45就是0.449近似而来的,表示第一组实验选择的硬币是A的概率为0.45。然后,0.449 * 5H = 2.2H ,0.449 * 5T = 2.2T ,表示第一组实验选择A硬币且正面朝上次数和反面朝上次数的期望值都是2.2,其他的值依次类推。最后,求出θ_A^{(1)}θA(1)=0.71,θ_B^{(1)}θB(1)=0.58。重复上述过程,不断迭代,直到算法收敛到一定精度为止。
这篇博客对EM算法的推导非常详细,链接如下:
https://blog.csdn.net/zhihua_oba/article/details/73776553
Python实现
#coding=utf-8
from numpy import *
from scipy import stats
import time
start = time.perf_counter()
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 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 iterations:最大迭代次数
:return:局部最优的模型参数
"""
iteration = 0;
while iteration < iterations:
new_prior = em_single(prior,observations)
delta_change = abs(prior[0]-new_prior[0])
if delta_change < tol:
break
else:
prior = new_prior
iteration +=1
return [new_prior,iteration]
#硬币投掷结果
observations = 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]])
print (em(observations,[0.6,0.5]))
end = time.perf_counter()
print('Running time: %f seconds'%(end-start))
来源:https://www.pianshen.com/article/7293314043/
猜你喜欢
- --利用T-SQL语句,实现数据库的备份与还原的功能 ----体现了SQL Server中的四个知识点: ----1. 获取SQL Serv
- 今天在玩 google earth 4.0b,发现 Print Screen 下来的图片很大,如果直接放在网页上,因为尺寸太大又不合适,又不
- 一、继承的格式类的继承格式如下,括号中的为父类名。class 类名(父类名): 代码二、单继承下面我们让Man继承自Master类,并调用继
- 学习目的: 学习ADO.NET用法,并如何用DataRearder读取数据 今天练习数据库的最基本用法,如何打开数据库。首先在网站设置文件w
- 本文实例讲述了Django框架静态文件处理、中间件、上传文件操作。分享给大家供大家参考,具体如下:Django静态文件处理、中间件、上传文件
- rem ---表单提示函数 Being-----------------------------CODE Copy ... Fu
- 前言:泛型是静态类型语言的基本特征,允许将类型作为参数传递给另一个类型、函数、或者其他结构。TypeScript 支持泛型作为将类型安全引入
- 代码如下:Set Catalog_object= Server.CreateObject("ADO
- 武器档案名称:firebug最新版本:1.7用途:前端调试器必备指数:使用难度:firebug是前端最具盛名的调试器,功能非常强悍。fire
- Python实现对变位词的判断,供大家参考,具体内容如下什么是变位词呢?即两个单词都是由相同的字母组成,而各自的字母顺序不同,譬如pytho
- 公布到网页上的Email经常会被一些工具自动提取,一些非法用户就会利用所提取的Email大肆发送垃圾邮件。这些工具大多都是查找链接中“mai
- 利用系统crontab来定时执行备份文件,按日期对备份结果进行保存,达到备份的目的。1、创建保存备份文件的路径/mysqldata#mkdi
- 目录:分析和设计组件编码实现和算法用 Ant 构建组件测试 JavaScript 组件我们走到哪儿了?前两期思考了太多东西,你是否已有倦意?
- 本文实例为大家分享了Python/C++实现字符串逆序的具体代码,供大家参考,具体内容如下题目描述:将字符串逆序输出Python实现一:借助
- 简述和GNU一样,YAML是一个递归着说“不”的名字。不同的是,GNU对UNIX说不,YAML说不的对象是XML。YAML不是XML。为什么
- 0 引言上周被一则新闻震惊到了,《2454万元大奖无人认领!福彩史上第二大弃奖在广东中山产生 》,在2019年5月2日开奖的双色球中,广东中
- 执行时间方法1import datetimestarttime = datetime.datetime.now()#long running
- 长话短说,今天介绍实现此功能的一个方法,需要了解的朋友可以参考下:一、JS 重载页面,本地刷新,返回上一页 代码如下:<a href=
- Div+CSS+JS组和能够实现很多好看的特殊的效果,这里推荐一款可刷新的下拉菜单:下面是js代码部分:<script type=te
- Python字典是另一种可变容器模型,且可存储任意类型对象,如字符串、数字、元组等其他容器模型。一、创建字典字典由键和对应值成对组成。字典也