注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

激光遥感

奋斗,探求,不达目的,誓不罢休!

 
 
 

日志

 
 
关于我

I was a university teacher, main interest areas:1、Spatio-temporal data analysis 2、Machine Learning3、 Pattern Recognition

网易考拉推荐

Metropolis Hasting算法  

2013-10-12 17:41:58|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

Metropolis Hasting Algorithm:

MH算法也是一种基于模拟的MCMC技术,一个很重要的应用是从给定的概率分布中抽样。主要原理是构造了一个精妙的Markov链,使得该链的稳态 是你给定的概率密度。它的好处,不用多说,自然是可以对付数学形式复杂的概率密度。有人说,单维的MH算法配上Gibbs Sampler几乎是“无敌”了。

今天试验的过程中发现,MH算法想用好也还不简单,里面的转移参数设定就不是很好弄。即使用最简单的高斯漂移项,方差的确定也是个头疼的问题,需要不同问题不同对待,多试验几次。当然你也可以始终选择“理想”参数。

还是拿上次的混合高斯分布来做模拟,模拟次数为500000次的时候,概率分布逼近的程度如下图。虽然几个明显的"峰"已经出来了,但是数值上还是 有很大差异的。估计是我的漂移方差没有选好。感觉还是inverse sampling好用,迭代次数不用很多,就可以达到相当的逼近程度。

试了一下MH算法,

R Code:

p=function(x,u1,sig1,u2,sig2){
(1/3)*(1/(sqrt(2*pi)*15)*exp(-0.5*(x-70)^2/15^2)+1/(sqrt(2*pi)*11)*exp(-0.5*(x+80)^2/11^2)+1/(sqrt(2*pi)*sig1)*exp(-0.5*(x-u1)^2/sig1^2)+1/(sqrt(2*pi)*sig2)*exp(-0.5*(x-u2)^2/sig2^2))
}


MH=function(x0,n){
x=NULL
x[1] = x0
for (i in 1:n){
x_can= x[i]+rnorm(1,0,3.25)
d= p(x_can,10,30,-10,10)/p(x[i],10,30,-10,10)
alpha= min(1,d)
u=runif(1,0,1)
if (u<alpha){
x[i+1]=x_can}
else{
x[i+1]=x[i]
}
if (round(i/100)==i/100) print(i)
}
x
}
z=MH(10,99999)
z=z[-10000]
a=seq(-100,100,0.2)

plot(density(z),col=1,main='Estimated Density',ylim=c(0,0.02),lty=1)
points(a, p(a,10,30,-10,10),pch='.',col=2,lty=2)
legend(60,0.02,c("True","Sim (MH)"),col=c(1,2),lty=c(1,2))

  评论这张
 
阅读(965)| 评论(0)
推荐

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017