前面已经大致的叙述了MCMC方法。今天来分享一下R中的一个实现MCMC算法的包mcmc。
mcmc包的一个核心函数就是metrop,其调用格式为:
metrop(obj, initial, nbatch, blen = 1, nspac = 1, scale = 1, outfun, debug = FALSE, ...)参数介绍:
obj:这里可以是一个metropo生成的对象,也可以是后验分布的对数似然函数
initial:初始值
nbatch:需要进行迭代的次数
blen:需要每批模拟的长度
这里metropo函数采用的模拟办法是正态随机游走,如果你要用其他办法去做,至少mcmc包没法帮你实现。
我们来看一个生成正态分布的例子:
h<-function(x){log(dnorm(x))}
out <- metrop(h, 0, 5000)
我们得到的模拟结果在out对象的batch中,除此之外,out还还包含了如下对象:
> names(out)
[1] "accept" "batch" "initial" "final"
[5] "initial.seed" "final.seed" "time" "lud"
[9] "nbatch" "blen" "nspac" "scale"
[13] "debug"
其中accept表示转移接受的概率,我们在之前的讨论就说过,一般accept要控制在40%左右为最佳,太大或太小对数据自相关性影响都比较大,过大或过小可以通过调整metropo中的scale参数得到改善。lud显示我们设置的对数似然。batch储存模拟结果。