原来的帖子见：http://www.madio.net/home.php?mod=space&uid=88948&do=blog&id=11156

 

function [mcs,efp]=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函数：

[y,rtp]=mhs(Normalpdf,50,1,1000)

初值为50，方差为1，长度为1000

运行后得到的链存在y中，可以做出y的路径图直方图：

plot(y)

hist(y,50)

结果如下：

可以看出，大概在520次转移时，链条开始收敛到标准正态分布

直方图：

这是整个链条的直方图，因为链条有约一半的数据还没有收敛到正态分布，所以不是很直观，如果把链条前面的550个数据“切掉”，只保留后面的450个数据明就好看多了，这也是为什么MCMC需要把“头”砍掉的原因

yhat=y([550:1000]);

hist(yhat,50)

这样结果就好多了。