R语言与Markov Chain Monte Carlo(MCMC)方法学习笔记(2)

   前面已经大致的叙述了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储存模拟结果。