贝叶斯集锦：R和JAGS的交互

转载自：http://site.douban.com/182577/widget/notes/10567181/note/295466672/

Markov chain Monte Carlo (MCMC)方法最早的实现是Linux下的BUGS，主要是用于Bayesian models涉及的统计计算（1989年），后来移植到Windows下发展成为WinBUGS，并终止了在Linux下的研发。它并不是开源的，于是芬兰的Helsinki大学搞了一个开源的OpenBUGS，法国人Martyn Plummer研发了个开源的JAGS。

JAGS，全称是Just another Gibbs sampler，是基于BUGS语言开发的利用MCMC来进行贝叶斯建模的软件包。它没有提供建模所用的GUI以及MCMC抽样的后处理，这些要在其它的程序软件上来处理，比如说利用R包（rjags）来调用JAGS并后处理MCMC的输出。JAGS相对于WinBUGS/OpenBUGS的主要优点在于平台的独立性，可以应用于各种操作系统，而WinBUGS/OpenBUGS只能应用于windows系统；JAGS也可以在64-bit平台上以64-bit应用来进行编译。

JAGS和R的交互非常好，下面我们使用rjags包来实现R对JAGS的调用。运行一个JAGS模型是指在参数的后验分布中生成抽样，需要这样5个步骤：

定义模型
初始化
编译
适应
监测
后续的MCMC收敛诊断、模型评价等等工作是要由R来完成的。当然，在使用rjags之前，要保证JAGS已经安装在你的电脑上。（JAGS下载：http://sourceforge.net/projects/mcmc-jags/files/）

下面采用对一个简单线性回归方程的仿真，来介绍用rjags包调用JAGS进行贝叶斯建模的过程。

对一个线性回归方程：
yi=α+βxi+ϵi

其中α，β是回归系数，ϵ是误差项，且ϵ∼N(0,σ2)。在线性回归中解释变量看做常量，而响应变量y视为随机变量。且y∼N(α+βx,σ2)。

- 用JAGS定义模型：

cat("model{for( i in 1:N ){y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + beta*x[i]}
alpha ~ dnorm(0, 1.0E-6)
beta ~ dnorm(0, 1.0E-6)
sigma ~ dunif(0,100)
tau <- 1/pow(sigma,2)}",file="F:/R/reg.jag")

JAGS的语法和R很相似。但对正态分布的密度dnorm来说，R的散布参数取成标准差，JAGS的散布参数取为方差的倒数，称为精度（precision），在上面的程序里我们把这个参数记为tau。建立模型之后，把它存入名为example1.jag的文件中。
三个参数$\alpha$，$\beta$，$\sigma$按无信息先验确定分布。$\alpha$，$\beta$的先验分布取成具有很大方差的正态分布，而$\sigma$取成一个较大区间中的均匀分布。

- 初始化

library(coda)
library(rjags)
reg.dat <- list( x=x, y=y, N=1 )
reg.ini <- list( list( alpha=0.05, beta=1.0, sigma=0.9 ),
list( alpha=0.04, beta=1.1, sigma=1.0 ),
list( alpha=0.06, beta=0.9, sigma=1.1 ) )
reg.par <- c("alpha","beta","sigma" )

数据reg.dat是一个列表，用来存放迭代中产生的x,y和N值。我们打算使用3条MCMC链来进行抽样（见下一段程序中的n.chains = 3），因此对每条链赋一个初始值，得到一个列表的列表reg.ini（如果赋成相同的初始值，那就是一个列表）。reg.par 里存放的是模型参数。

- 适应和编译

reg.mod <- jags.model( file = "F:/R/reg.jag",
data = reg.dat,
n.chains = 3,
inits = reg.ini,
n.adapt = 2000 )

###Compiling model graph
###Resolving undeclared variables
###Allocating nodes
###Graph Size: 3113

###Initializing model

- 监测和图形诊断

通过调用jags.model函数，写在m2.jag中的程序就被JAGS编译了，对未定义的参数和变量进行检查，做好进行后验抽样的准备。有时候模型中参数很多，需要指定进行抽样的参数，这个过程叫做检测（monitoring）。

给定一个jags.model对象，可以通过调用coda.samples函数来对抽样个数进行收集存为coda对象。然后利用plot函数绘出由coda包所提供的图形诊断。

reg.res <- coda.samples( reg.mod,
var = reg.par,
n.iter = 10000,
thin = 1)
plot(reg.res)

图形显示的是2000次MCMC迭代的后验抽样。

- burn-in和后验概括

因为MCMC抽样过程中，马尔科夫链是随n增大而收敛。在使用MCMC抽样的结果时，要有合适一个burn-in的过程，使得burn-in之后的链的样本是收敛的。
取burn-in值为1000，对之后的样本采用summary函数来得到后验抽样的统计量。

burn.in <- 5000
a<-window(reg.res, start = burn.in)
summary(window(reg.res, start = burn.in))

###Iterations = 5000:12000
###Thinning interval = 1
###Number of chains = 3
###Sample size per chain = 7001

###1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean SD Naive SE Time-series SE
alpha -1.712 57.96 0.3999 0.3986
beta 3.824 1001.52 6.9107 6.9511
sigma 50.203 28.85 0.1991 0.4586

###2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%
alpha -130.35 -27.75 -1.308 24.41 124.24
beta -1955.49 -683.66 12.810 682.41 1962.03
sigma 2.47 25.25 50.210 75.24 97.62

贝叶斯集锦：R和JAGS的交互相关推荐

贝叶斯集锦：MCMCpack包
转载自:http://site.douban.com/182577/widget/notes/10567181/note/280112466/ ####贝叶斯集锦这个系列目的是想收集一些使用R的贝叶斯 ...
贝叶斯集锦：贝叶斯统计基础
转载自:http://site.douban.com/182577/widget/notes/10567181/note/294041203/ 1.从贝叶斯定理到贝叶斯统计推断 (1)贝叶斯统计简史 ...
贝叶斯集锦：从MC、MC到MCMC
转载自: #####一份草稿贝叶斯计算基础一.从MC.MC到MCMC 斯坦福统计学教授Persi Diaconis是一位传奇式的人物.Diaconis14岁就成了一名魔术师,为了看懂数学家Fell ...
贝叶斯集锦：贝叶斯派和频率派的一个例子
转载自:http://site.douban.com/182577/widget/notes/10567181/note/278503359/ 这个例子的主要目的在于探讨贝叶斯派和频率派适用的具体情境 ...
机械学习与R语言--Naive Bayes 朴素贝叶斯在R语言中的实现
为什么天气预报说70%概率下雨?为什么垃圾短信垃圾邮件被自动归类?这一切的基础算法便是朴素贝叶斯理论(算法有很多,这仅是其中之一). 1.由贝叶斯理论到朴素贝叶斯(naive bayes) 理论的基础 ...
R语言贝叶斯广义线性混合（多层次/水平/嵌套）模型GLMM、逻辑回归分析教育留级影响因素数据...
全文下载链接:http://tecdat.cn/?p=24203 本教程使用R介绍了具有非信息先验的贝叶斯 GLM(广义线性模型) (点击文末"阅读原文"获取完整代码数据). 当前 ...
R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型
全文下载链接:http://tecdat.cn/?p=11974 R2WinBUGS软件包提供了从R调用WinBUGS的便捷功能.它自动以WinBUGS可读的格式写入数据和脚本,以进行批处理(自1.4 ...
概率图模型基于R语言(一)贝叶斯模型
概率图模型基于R语言[一]贝叶斯模型条件概率贝叶斯模型 R 参考书:<概率图模型基于R语言> 记录一些代码和自己的想法和目前版本代码的写法(书中有一些无法用了) 随时更新条件概率熟 ...
【译文】利用STAN做贝叶斯回归分析：Part 1 正态回归
[译文]利用STAN做贝叶斯回归分析:Part 1 正态回归作者 Lionel Hertzog 本文将介绍如何在R中做贝叶斯回归分析,你能在文末的参考文献中找到相关主题的更多信息. 贝叶斯回归贝叶 ...

贝叶斯集锦：R和JAGS的交互

贝叶斯集锦：R和JAGS的交互相关推荐

最新文章

热门文章