马尔科夫链蒙特卡洛（二）

​来，做一道书后习题：用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分布很像吧