QQ登录

只需要一步,快速开始

 注册地址  找回密码

tag 标签: MCMC

相关日志

分享 MCMC的MH抽样MATLAB程序说明
liwenhui 2013-5-30 15:06
MCMC的MH抽样MATLAB程序说明
原来的帖子见: http://www.madio.net/home.php?mod=spaceuid=88948do=blogid=11156 function =mhs(f,p0,sigma,rt) %mcs 为得到的马尔可夫链,它是一个列向量的形式给出 %efp 为程序中迭代的效率指标,它等于得到的链的长度比上总的迭代次数 %f 为目标函数,它是需要通过 MCMC 生成的目标概率密度函数,可以从 m 问文件定义,也可以用匿名函数定义 %p0 为马尔可夫链的初始值,理论上讲它可以为任何值,但是为了能让链条尽快收敛于平稳分布(即上面的 f ),它最好取接近平稳分布的均值附近。 P0 将会在初次迭代中充当转移核的均值,并依次不断被新的值替代,由此构造了一个条件分布函数 %sigma 转移核的方差,方差大的话收敛会快一点,但可能不稳定,方差小的话收敛会比较满,链条长度最好取长一些 %rt 指定目标链条的长度 k1=1;k2=0; gser=zeros(rt,1);gser(1)=p0; while (k1 rt); p1=p0+sigma*randn; q1=feval(f,p0); q0=feval(f,p1); r=q0/q1; alpha=min(1,r); urand=rand(1); if urand alpha; gser(k1+1)=p1;p0=p1;k1=k1+1; end k2=k2+1; end mcs=gser;efp=rt/k2; -------------------------------------------------------------------------------------------------------------------------------------- 将以上程序存为 m 文件,命名为 mhs.m ,调用时先定义目标平稳分布的密度函数 f ,然后调用 mhs ,指定 p0 ,sigma,rt. 比如要生成一个符合标准正态分布的马尔可夫链(实践中一般不会通过 MCMC 来生成正态分布,因为正态分布有其它更好的方法,此处只是举个例子) 首先定义正态分布函数 : Normalpdf=@(x)((1./sqrt(2.*pi)).*exp(-0.5.*x.^2)) 然后调用 mhs 函数: =mhs(Normalpdf,50,1,1000) 初值为 50 ,方差为 1 ,长度为 1000 运行后得到的链存在 y 中,可以做出 y 的路径图直方图: plot(y) hist(y,50) 结果如下: 可以看出,大概在 520 次转移时,链条开始收敛到标准正态分布 直方图: 这是整个链条的直方图,因为链条有约一半的数据还没有收敛到正态分布,所以不是很直观,如果把链条前面的 550 个数据“切掉”,只保留后面的 450 个数据明就好看多了,这也是为什么 MCMC 需要把“头”砍掉的原因 yhat=y( ); hist(yhat,50) 这样结果就好多了。
个人分类: 在人间|2591 次阅读|0 个评论
qq
收缩
  • 电话咨询

  • 04714969085

关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

手机版|Archiver| |繁體中文 手机客户端  

蒙公网安备 15010502000194号

Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

GMT+8, 2025-7-20 22:28 , Processed in 0.249167 second(s), 24 queries .

回顶部