2011年1月11日星期二

马尔科夫链蒙特卡洛(二)

马尔科夫链蒙特卡洛(二)

来,做一道书后习题:用MCMC生成服从Poisson分布的随机数。

习题嘛,瞎折腾,本来R就有现成的函数:rpois
​既然MCMC,首先要构造一个转移规则。考虑射线上的随机游动。从0出发,必然走到1,从k出发(k>0),一半的概率向左右一半的概率向右走。

计算接受概率:

其中

其中k>1。边界0和1处需要稍微特殊处理一下。注意到,概率密度函数求比值之后,我们再也不用算指数和阶乘了。剩下的事情,就是不停地走,走每一步时只需要回答:

1)往哪边走?---用随机游动的trial distribution:一半的概率向左右一半的概率向右走

2)真的要走到那里去吗?---使用上面算的接受概率(MCMC的核心)

下面是R源代码


p=function(s){
    runif(1) <= s
}
#pi_rate(k,lambda)=pi(k,lambda)/pi(k+1,lambda)
pi_rate=function(k,lambda){
    (k+1)/lambda
}
accept_walk_left_from=function(k,lambda){
    if (k==1)
        min(1, 2 / lambda)
    else
        min(1, pi_rate(k - 1, lambda))
}
accept_walk_right_from=function(k, lambda){
    if (k==0)
        min(1, lambda / 2)
    else
        min(1, 1 / pi_rate(k, lambda))
}
transition_from=function(k, lambda){
    if (k==0){
        if (p(accept_walk_right_from(k, lambda)))
            1
        else
            0
    }
    else{
        if (p(0.5)){
            if (p(accept_walk_right_from(k, lambda)))
                k + 1
            else
                k
        }else{
            if (p(accept_walk_left_from(k, lambda)))
                k - 1
            else
                k
        }
    }
}
mcmc_poisson=function(n,lambda){
    f=0
    for(i in 1:n){
        f[i+1]=transition_from(f[i], lambda)
    }
    f
}

samples=tail(mcmc_poisson(20000, 5), 10000)
hist(samples)


下面是一次实验结果的直方图,看起来和Poisson分布很像吧




没有评论: