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)