一、EM算法

EM算法是一种在观测到数据后,用迭代法估计未知参数的方法。可以证明EM算法得到的序列是稳定单调递增的。这种算法对于截尾数据或参数中有一些我们不感兴趣的参数时特别有效。

EM算法的步骤为:

E-step(求期望):在给定y及theta=theta(i)的条件下,求关于完全数据对数似然关于潜在变量z的期望

M-step(求极值):求上述期望关于theta的最大值theta(i+1)

重复以上两步,直至收敛即可得到theta的MLE。

从上面的算法我们可以看到对于一个参数的情况,EM仅仅只是求解MLE的一个迭代算法。M-step做得就是optimize函数做得事情。对于EM算法,我们也没有现成的求解函数(这个是自然的),我们一样可以通过人机交互的办法处理。

先举一个一元的例子:

设一次实验可能有4个结果,发生概率分别为0.5+theta/4, 0.25-theta/4 ,0.25-theta/4 ,theta/4.其中theta在0,1之间。现进行了197次实验,结果发生的次数分别为:125,18,20,34,求theta的MLE。

计算出theta(i+1)=(195theta(i)+68)/(197theta(i)+144)

为什么是这个结果,请翻阅王兆军《数理统计讲义》p43-p44

我们用简单的循环就可以解决这个问题,程序及结果如下:

>fun<-function(error=1e-7){

+theta<-0.5

+k<-1

+while(T){

+k<-k+1

+theta[k]<-(159*theta[k-1]+68)/(197*theta[k-1]+144)

+if(abs(theta[k]-theta[k-1])<error) break

+}

+list(theta<-theta[k],iter<-k)

+}

>fun()

[[1]]

[1]0.6268215

[[2]]

[1]9

我们再看一个二元的简单例子:

幼儿园里老师给a,b,c,d四个小朋友发糖吃,但老师有点偏心,不同小朋友得到糖的概率不同,p(a)=0.5,p(b)=miu, p(c)=2*miu, p(d)=0.5-3*miu 如果确定了参数miu,概率分布就知道了。我们可以通过观察样本数据来推测参数知道c和d二人得到的糖果数,也知道a与b二人的糖果数之和为h,如何来估计出参数miu呢?前面我们知道了,如果观察到a,b,c,d就可以用ML估计出miu。反之,如果miu已知,根据概率期望 a/b=0.5/miu,又有a+b=h。由两个式子可得到  a=0.5*h/(0.5+miu)和b=miu*h/(0.5+miu)。

># 已知条件

>

>h = 20

>c = 10

>d = 10

>

># 随机初始两个未知量

>miu = runif(1,0,1/6)

>b = round(runif(1,1,20))

>

>iter = 1

>nonstop=TRUE

>while (nonstop) {

+     # E步骤,根据假设的miu来算b

+     b = c(b,miu[iter]*h/(0.5+miu[iter]))

+     print(b)

+     # M步骤,根据上面算出的b再来计算miu

+     miu = c(miu,(b[iter+1] +c)/(6*(b[iter+1]+c+d)))

+     print(miu)

+     # 记录循环次数

+     iter = iter + 1

+     # 如果前后两次的计算结果差距很小则退出

+     nonstop =((miu[iter]-miu[iter-1])>10^(-10))

+}

[1]3.000000 4.450531

[1]0.14310878 0.09850182

>print(cbind(miu,b))

miu        b

[1,]0.14310878 3.000000

[2,]0.09850182 4.450531

关于EM算法,及后续的发展GME的理论你可以在多数数理统计书上找到相关结论,也可以用类似办法编写函数处理它。

二、       自助法(bootstrap)

Bootstrap法是以原始数据为基础的模拟抽样统计推断法,可用于研究一组数据的某统计量的分布特征,特别适用于那些难以用常规方法导出对参数的区间估计、假设检验等问题。“Bootstrap”的基本思想是:在原始数据的围内作有放回的再抽样,样本含量仍为n,原始数据中每个观察单位每次被抽到的概率相等,为1,…,n,所得样本称为bootstrap样本。于是可得到参数Η的一个估计值Η(b),这样重复若干次,记为B。设B=1000,就得到该参数的1000个估计值,则参数Η的标准误的bootstrap估计。简而言之就是:既然样本是抽出来的,那我何不从样本中再抽样。

我们知道,如果分布函数F是已知的。在理论上就能够计算出参数的估计量的均方误差.若分布函数f未知,由格里文科-康特利定理知,当M充分大时,经验分布函数以概率1一致收敛到F。

我们举一例:利用bootstrap法估计标准正态分布随机变量的期望theta=E(X)

>gauss<-rnorm(100,2,6)

>boot<-0

>for(i in 1:1000){

+boot[i]=mean(sample(gauss,replace=T))

+}

>summary(boot)

Min. 1st Qu. Median    Mean 3rd Qu.    Max.

0.3345 1.9540  2.3350  2.3230 2.7020  4.2330

>summary(gauss)

Min. 1st Qu. Median    Mean 3rd Qu.    Max.

-13.380  -2.238  2.570   2.296   6.861 16.230

>sd(boot)

[1]0.599087

>sd(gauss)/sqrt(100)

[1]0.5906275

结果分析:

方法 标准差 均值
理论值 0.6 2
矩估计 0.5906275 2.570
bootstrap算法 0.599087 2.3230

需要指出的是bootstrap法不是为了提高估计量的精度.而是一般用来对估计量的方差进行估计。

关于Bootstrap方法,R中有相应的扩展包,你可以通过??bootstrap来查看有关内容

R语言与点估计学习笔记(EM算法与Bootstrap法)相关推荐

  1. R语言与点估计学习笔记(矩估计与MLE)

    众所周知,R语言是个不错的统计软件.今天分享一下利用R语言做点估计的内容.主要有:矩估计.极大似然估计.EM算法.最小二乘估计.刀切法(Jackknife).自助法(Bootstrap)的相关内容. ...

  2. R语言与点估计学习笔记(刀切法与最小二乘估计)

    一.       刀切法(jackknife) 刀切法的提出,是基于点估计准则无偏性.刀切法的作用就是不断地压缩偏差.但需要指出的是缩小偏差并不是一个好的办法,因为偏差趋于0时,均方误差会变得十分大. ...

  3. R语言与机器学习学习笔记(分类算法)

    转载自:http://www.itongji.cn/article/0P534092014.html 人工神经网络(ANN),简称神经网络,是一种模仿生物神经网络的结构和功能的数学模型或计算模型.神经 ...

  4. 【统计学习方法】学习笔记——EM算法及其推广

    统计学习方法学习笔记--EM算法及其推广 1. EM算法的引入 1.1 EM算法 1.2 EM算法的导出 1.3 EM算法在非监督学习中的应用 2. EM算法的收敛性 3. EM算法在高斯混合模型学习 ...

  5. 语言 提取列名_学习健明老师发布的R语言练习题的学习笔记(二)

    学习者:骆栢维 题目来源:生信基石之R语言 中级10 个题目:http://www.bio-info-trainee.com/3750.html 备注:本文为笔者学习健明老师GitHub答案代码的学习 ...

  6. 【转载】R语言与数据挖掘学习笔记

    (1):数据挖掘相关包的介绍 今天发现一个很不错的博客(http://www.RDataMining.com),博主致力于研究R语言在数据挖掘方面的应用,正好近期很想系统的学习一下R语言和数据挖掘的整 ...

  7. R plot图片背景设置为透明_学习健明老师发布的R语言练习题的学习笔记(一)...

    学习者:骆栢维 题目来源:生信基石之R语言 初级10 个题目:http://www.bio-info-trainee.com/3793.html 备注:本文为笔者学习健明老师GitHub答案代码的学习 ...

  8. R语言dplyr包学习笔记(吐血整理宇宙无敌详细版)

    出处:AI入门学习 dplyr包主要用于数据清洗和整理,主要功能有:行选择.列选择.统计汇总.窗口函数.数据框交集等是非常高效.友好的数据处理包,学清楚了,基本上数据能随意玩弄,对的,随意玩弄,简直大 ...

  9. R语言基础知识-学习笔记汇总

    B站课程:生信必备技巧之R语言基础教程全集的代码笔记 1.R语言包安装 rm(list = ls()) # 设置镜像: options()$repos options()$BioC_mirror #o ...

最新文章

  1. html5 点击事件委托,jquery事件委托
  2. linux ext4 img解包打包教程,解打包.img.ext4(转)
  3. 51nod1228 序列求和(伯努利数)
  4. C#垃圾回收学习总结
  5. 简单多进程任务处理程序
  6. java中xpath_java-xpath学习
  7. 【Python基础】Pandas批量合并文件脚本,多个同名sheet也适用
  8. AspNetPager分页控件的运用 【转】-有用
  9. 信息化工程监理规范_房建工程监理资料管理存在的问题及应对措施
  10. windows无法格式化u盘_windows无法完成格式化怎么办
  11. mysql 循环 索引值,mysql:循环遍历表和alter table添加索引
  12. excel计算机快捷键大全,excel表格使用技巧快捷键大全
  13. VS2010 上手案例---hello word
  14. 打印机服务器没有响应 请检查设置,打印机服务无法启动的解决办法
  15. 30天自制操作系统——第一天到第二天
  16. CS:APP Archlab(未完待续)
  17. 数字IC面试高频考点之跨时钟域信号处理
  18. WPS的Excel做一个下拉选择功能
  19. 微信小程序——增删改
  20. 加入域时出现“不能访问网络位置”错误信息

热门文章

  1. CS品牌SD NAND在点读机产品中的案例分享
  2. Nginx学习笔记5--(极客时间-陶辉)
  3. 【日语】标日初级上册单词(1-4)1
  4. Zynq移植USB触摸屏
  5. 骚扰电话短信、诈骗电话拦截
  6. 公司来了一个女程序员,新鲜
  7. 上海将设大数据交易中心 市场化民营化运作
  8. Perfdump 工具
  9. 被傅里叶变换,劝退的著名学科导师
  10. 自我管理方法办法概要