用MATLAB编写MCMC算法产生一条平稳分布为指定分布的马尔可夫链：

function [mcs,efp]=mhs(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;

f可以用M文件定义，也可以用匿名函数定义。