R语言:R2OpenBUGS
用这个包调用BUGS model,分别用表格和图形概述inference和convergence,保存估计的结果
as.bugs.array 转换成bugs object
函数把马尔科夫链估计结果(不是来自于BUGS),转成BUGS object,主要用来plot.bugs 展示结果。
as.bugs.array(sims.array, model.file=NULL, program=NULL, DIC=FALSE, DICOutput=NULL, n.iter=NULL, n.burnin=0, n.thin=1)
sims.array :3维数组 n.keep, n.chains和combined parameter vector的长度
model.file : OpenBUGS编写的.odc 模型文件
DIC : 是否计算DIC曲线
DICOutput : DIC值
n.iter :生成sims.array 每条chain 迭代数
n.burnin :丢弃的迭代次数
n.thin :thinning rate
attach.all 添加数据到搜索路径
The database is attached/detached to the search path,While attach.all attaches all elements of an object x to a database called name, attach.bugs attaches all elements of x$sims.list to the database bugs.sims itself making use of attach.all.
attach.all(x, overwrite = NA, name = “attach.all”) attach.bugs(x, overwrite = NA) detach.all(name = “attach.all”) detach.bugs()
x : bugs 对象
overwrite :TRUE 删除全局环境中被掩盖的数据, NA 询问,FALSE
name : 环境name
bugs 最重要,用R运行bugs
自动输入值,启动bugs,保存结果
bugs(data, inits, parameters.to.save, n.iter, model.file=“model.txt”, n.chains=3, n.burnin=floor(n.iter / 2), n.thin=1, saveExec=FALSE,restart=FALSE, debug=FALSE, DIC=TRUE, digits=5, codaPkg=FALSE, OpenBUGS.pgm=NULL, working.directory=NULL, clearWD=FALSE, useWINE=FALSE, WINE=NULL, newWINE=TRUE, WINEPATH=NULL, bugs.seed=1, summary.only=FALSE, save.history=(.Platform$OS.type == “windows” | useWINE==TRUE), over.relax = FALSE)
data :模型中使用的数据
inits :n chain 的元素列表,每一个要素是一个模型初始值列表,或者一个生成初始值得function
parameters.to.save : 需要被记录的参数名向量
model.file : model 文件.txt
n.chains : 默认3条
n.iter :每条链的迭代次数,默认2000
n.thin : Thinning rate. 正整数,默认是1,
saveExec :使用basename(模型.file)保存OpenBUGS执行的重新启动映像。
restart :执行从上次执行的最后状态恢复,存储在工作目录中的.bug文件中。
debug : 默认FALSE,正在运行行时Openbugs 页面关闭
DIC :计算deviance,,pD,和DIC。
digits :有效小数位数
codaPkg :FALSE 返回 bugs对象,否则输出,用coda 包 read.bugs 读取,
OpenBUGS.pgm :通向OpenBUGS可执行程序的完整路径。
working.directory:OpenBUGS的输入和输出将存储在此目录中;
clearWD :是否这些文件的“data.txt”、“init(1:n.chains). txt”,“log.odc”、“codaIndex.txt”和“coda[1:nchain].txt”结束时删除。
bugs.seed :OpenBUGS随机种子,1-14整数
summary.only : TURE ,仅对非常快速的分析给出了一个参数概要
save.history : TURE,最后画出trace
# An example model file is given in:
model.file <- system.file(package="R2OpenBUGS", "model", "schools.txt")
# Let's take a look:
print(model.file)
[1] "C:/Users/Date/Documents/R/win-library/3.5/R2OpenBUGS/model/schools.txt"
file.show(model.file)
model { for (j in 1:J){ y[j] ~ dnorm (theta[j], tau.y[j]) theta[j] ~ dnorm (mu.theta, tau.theta) tau.y[j] <- pow(sigma.y[j], -2) } mu.theta ~ dnorm (0.0, 1.0E-6) tau.theta <- pow(sigma.theta, -2) sigma.theta ~ dunif (0, 1000) }
data(schools)
schools
school estimate sd
1 A 28.39 14.9
2 B 7.94 10.2
3 C -2.75 16.3
4 D 6.82 11.0
5 E -0.64 9.4
6 F 0.63 11.4
7 G 18.01 10.4
8 H 12.16 17.6
J <- nrow(schools)
y <- schools$estimate
sigma.y <- schools$sd
data <- list ("J", "y", "sigma.y")
data
[[1]]
[1] "J"
[[2]]
[1] "y"
[[3]]
[1] "sigma.y"
inits <- function(){
list(theta=rnorm(J, 0, 100), mu.theta=rnorm(1, 0, 100),sigma.theta=runif(1, 0, 100))
}
parameters <- c("theta", "mu.theta", "sigma.theta")
schools.sim <- bugs(data, inits, parameters, model.file,
n.chains=3, n.iter=5000)
print(schools.sim)
Inference for Bugs model at "C:/Users/Date/Documents/R/win-library/3.5/R2OpenBUGS/model/schools.txt",
Current: 3 chains, each with 5000 iterations (first 2500 discarded)
Cumulative: n.sims = 7500 iterations saved
mean sd 2.5% 25% 50% 75% 97.5%
theta[1] 11.2 8.9 -2.8 5.5 9.8 15.7 32.9
theta[2] 7.5 6.5 -4.9 3.4 7.4 11.5 20.8
theta[3] 5.8 8.0 -12.3 1.3 6.4 10.5 21.0
theta[4] 7.1 6.7 -6.3 3.0 7.1 11.2 20.9
theta[5] 4.9 6.3 -8.8 1.0 5.5 9.1 16.4
theta[6] 5.8 6.8 -9.2 1.7 6.1 10.1 18.3
theta[7] 10.4 7.2 -2.1 5.6 9.7 14.6 26.3
theta[8] 8.1 8.0 -7.1 3.4 7.8 12.4 25.7
mu.theta 7.6 5.4 -2.7 4.1 7.6 11.0 18.6
sigma.theta 6.7 5.9 0.2 2.3 5.3 9.4 21.8
deviance 60.5 2.3 57.0 59.2 60.1 61.6 66.2
Rhat n.eff
theta[1] 1 410
theta[2] 1 420
theta[3] 1 1200
theta[4] 1 610
theta[5] 1 680
theta[6] 1 490
theta[7] 1 300
theta[8] 1 420
mu.theta 1 310
sigma.theta 1 540
deviance 1 1600
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.9 and DIC = 63.4
DIC is an estimate of expected predictive error (lower deviance is better).
plot(schools.sim)
bugs.data 生成输入文件
bugs.data(data, dir = getwd(), digits = 5, data.file = “data.txt”)
bugs.inits 生成初始值文件
bugs.inits(inits, n.chains, digits, inits.files = paste(“inits”, 1:n.chains, “.txt”, sep = “”))
bugs.log 读取log文件(summary statistics and DIC information)
bugs.log(file)
plot.bugs 画bugs对象
plot(x, display.parallel = FALSE, …)
display.parallel :在摘要图的两部分中显示平行的间隔
print.bugs 输出bugs对象
print(x, digits.summary = 1, …)
digits.summary:四舍五入的位数
read.bugs
读Markov链蒙特卡罗输出的CODA格式。并返回一个类mcmc.list对象。使用coda包进行进一步的输出分析列表。
read.bugs(codafiles, …)
validateInstallOpenBUGS 比较R和openbugs软件运行结果
write.model 转化R function创建模型文件
schoolsmodel <- function(){
for (j in 1:J){
y[j] ~ dnorm (theta[j], tau.y[j])
theta[j] ~ dnorm (mu.theta, tau.theta)
tau.y[j] <- pow(sigma.y[j], -2)
}
mu.theta ~ dnorm (0.0, 1.0E-6)
tau.theta <- pow(sigma.theta, -2)
sigma.theta ~ dunif (0, 1000)
}
## some temporary filename:
filename <- file.path(tempdir(), "schoolsmodel.txt")
## write model file:
write.model(schoolsmodel, filename)
## and let's take a look:
file.show(filename)