R2WinBUGSとR2OpenBUGSとでは、bugs()の返り値がちがう [統計]
R2OpenBUGSのbugs()関数は、n.thinによらず、すべてのマルコフ鎖の値を保持している。
コード
結果
コード
library(R2WinBUGS)
model.file <- system.file(package = "R2WinBUGS", "model", "schools.txt")
data(schools)
J <- nrow(schools)
y <- schools$estimate
sigma.y <- schools$sd
data <- list ("J", "y", "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")
WINE <- "/Applications/Darwine/Wine.bundle/Contents/bin/wine"
WINEPATH <- "/Applications/Darwine/Wine.bundle/Contents/bin/winepath"
post.winbugs <- bugs(data, inits, parameters, model.file,
n.chains = 3, n.iter = 11000, n.burnin = 1000, n.thin = 10,
working.directory = NULL, clearWD = TRUE,
useWINE = TRUE, newWINE = TRUE,
WINE = WINE, WINEPATH = WINEPATH, debug = FALSE)
print(post.winbugs)
str(post.winbugs)
#
library(R2OpenBUGS)
openbugs.version <- '321'
openbugs <- paste(Sys.getenv("HOME"),
"/.wine/drive_c/Program Files/OpenBUGS/OpenBUGS",
openbugs.version, "/OpenBUGS.exe", sep = "")
post.openbugs <- bugs(data, inits, parameters, model.file,
n.chains = 3, n.iter = 11000, n.burnin = 1000, n.thin = 10,
working.directory = NULL, clearWD = TRUE,
OpenBUGS.pgm = openbugs,
useWINE = TRUE, newWINE = TRUE,
WINE = WINE, WINEPATH = WINEPATH, debug = FALSE)
print(post.openbugs)
str(post.openbugs)
結果
> print(post.winbugs) Inference for Bugs model at "/Library/Frameworks/R.framework/Versions/2.13/Resources/library/R2WinBUGS/model/schools.txt", fit using WinBUGS, 3 chains, each with 11000 iterations (first 1000 discarded), n.thin = 10 n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff theta[1] 11.6 8.4 -1.8 6.2 10.3 15.8 31.8 1 500 theta[2] 8.1 6.3 -5.0 4.2 7.9 11.9 21.4 1 850 theta[3] 6.5 7.7 -10.6 2.5 7.0 11.2 20.6 1 970 theta[4] 7.8 6.5 -6.1 4.2 7.7 11.8 20.8 1 790 theta[5] 5.7 6.5 -9.3 2.1 6.2 9.9 17.4 1 1600 theta[6] 6.4 7.0 -9.2 2.4 6.8 10.8 19.4 1 3000 theta[7] 10.7 6.8 -0.8 6.4 10.0 14.5 26.1 1 690 theta[8] 8.7 7.6 -6.4 4.2 8.3 12.9 24.8 1 1900 mu.theta 8.4 5.1 -1.1 5.1 8.2 11.5 18.5 1 970 sigma.theta 6.5 5.7 0.2 2.4 5.3 8.9 21.2 1 120 deviance 60.5 2.2 57.1 59.2 60.1 61.4 66.0 1 3000 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.7 and DIC = 63.2 DIC is an estimate of expected predictive error (lower deviance is better). > str(post.winbugs) List of 24 $ n.chains : num 3 $ n.iter : num 11000 $ n.burnin : num 1000 $ n.thin : num 10 $ n.keep : num 1000 $ n.sims : num 3000 $ sims.array : num [1:1000, 1:3, 1:11] 13.8 4.2 4.43 10.83 -2.8 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : NULL .. ..$ : NULL .. ..$ : chr [1:11] "theta[1]" "theta[2]" "theta[3]" "theta[4]" ... $ sims.list :List of 4 ..$ theta : num [1:3000, 1:8] 12.93 5.04 13.13 13.38 4.56 ... ..$ mu.theta : num [1:3000] 10.8 12.1 12.7 14.9 6.5 ... ..$ sigma.theta: num [1:3000] 7.91 9.71 1.04 9.27 2.52 ... ..$ deviance : num [1:3000] 63.5 58.5 62 60.2 59.8 ... $ sims.matrix : num [1:3000, 1:11] 12.93 5.04 13.13 13.38 4.56 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:11] "theta[1]" "theta[2]" "theta[3]" "theta[4]" ... $ summary : num [1:11, 1:9] 11.58 8.08 6.55 7.75 5.74 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:11] "theta[1]" "theta[2]" "theta[3]" "theta[4]" ... .. ..$ : chr [1:9] "mean" "sd" "2.5%" "25%" ... $ mean :List of 4 ..$ theta : num [1:8(1d)] 11.58 8.08 6.55 7.75 5.74 ... ..$ mu.theta : num 8.38 ..$ sigma.theta: num 6.5 ..$ deviance : num 60.5 $ sd :List of 4 ..$ theta : num [1:8(1d)] 8.39 6.28 7.7 6.52 6.53 ... ..$ mu.theta : num 5.05 ..$ sigma.theta: num 5.69 ..$ deviance : num 2.21 $ median :List of 4 ..$ theta : num [1:8(1d)] 10.34 7.92 7 7.74 6.23 ... ..$ mu.theta : num 8.15 ..$ sigma.theta: num 5.28 ..$ deviance : num 60.1 $ root.short : chr [1:4] "theta" "mu.theta" "sigma.theta" "deviance" $ long.short :List of 4 ..$ : int [1:8] 4 5 6 7 8 9 10 11 ..$ : int 2 ..$ : int 3 ..$ : int 1 $ dimension.short: num [1:4] 1 0 0 0 $ indexes.short :List of 4 ..$ :List of 8 .. ..$ : num 1 .. ..$ : num 2 .. ..$ : num 3 .. ..$ : num 4 .. ..$ : num 5 .. ..$ : num 6 .. ..$ : num 7 .. ..$ : num 8 ..$ : NULL ..$ : NULL ..$ : NULL $ last.values :List of 3 ..$ :List of 3 .. ..$ theta : num [1:8] 12.42 7.26 9.18 19.57 11.64 ... .. ..$ mu.theta : num 11.5 .. ..$ sigma.theta: num 4.54 ..$ :List of 3 .. ..$ theta : num [1:8] 7.19 1.32 13.2 3.9 2.96 ... .. ..$ mu.theta : num 6.29 .. ..$ sigma.theta: num 6.34 ..$ :List of 3 .. ..$ theta : num [1:8] 7.34 6.88 7.55 7.28 8.48 ... .. ..$ mu.theta : num 7.22 .. ..$ sigma.theta: num 0.895 $ isDIC : logi TRUE $ DICbyR : logi FALSE $ pD : num 2.74 $ DIC : num 63.2 $ model.file : chr "/Library/Frameworks/R.framework/Versions/2.13/Resources/library/R2WinBUGS/model/schools.txt" $ program : chr "WinBUGS" - attr(*, "class")= chr "bugs" > print(post.openbugs) Inference for Bugs model at "/Library/Frameworks/R.framework/Versions/2.13/Resources/library/R2WinBUGS/model/schools.txt", Current: 3 chains, each with 11000 iterations (first 1000 discarded), n.thin = 10 Cumulative: n.sims = 30000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff theta[1] 11.8 8.4 -1.9 6.3 10.6 16.0 31.9 1 26000 theta[2] 8.1 6.4 -4.7 4.0 8.1 12.1 21.0 1 23000 theta[3] 6.3 7.9 -11.7 2.1 6.8 11.2 20.9 1 21000 theta[4] 7.8 6.6 -5.9 3.6 7.8 11.9 20.8 1 30000 theta[5] 5.5 6.5 -8.7 1.6 5.9 9.8 17.2 1 30000 theta[6] 6.2 6.9 -8.9 2.2 6.7 10.7 18.8 1 30000 theta[7] 10.8 6.9 -1.5 6.2 10.2 14.7 26.3 1 30000 theta[8] 8.7 7.9 -7.0 4.0 8.5 13.0 25.7 1 30000 mu.theta 8.1 5.3 -2.2 4.8 8.1 11.4 18.6 1 30000 sigma.theta 6.6 5.7 0.2 2.5 5.3 9.2 20.9 1 3900 deviance 60.4 2.2 57.0 59.1 60.0 61.4 66.0 1 7900 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 = 63.2 and DIC = 2.8 DIC is an estimate of expected predictive error (lower deviance is better). > str(post.openbugs) List of 23 $ n.chains : num 3 $ n.iter : num 11000 $ n.burnin : num 1000 $ n.thin : num 10 $ n.keep : num 10000 $ n.sims : num 30000 $ sims.array : num [1:10000, 1:3, 1:11] 10.82 8.92 5.37 13.24 5.3 ... ..- attr(*, "dimnames")=List of 3 .. ..$ : NULL .. ..$ : NULL .. ..$ : chr [1:11] "theta[1]" "theta[2]" "theta[3]" "theta[4]" ... $ sims.list :List of 4 ..$ theta : num [1:30000, 1:8] 5.65 12.91 9.5 6.81 5.63 ... ..$ mu.theta : num [1:30000] 7.616 7.288 7.762 8.155 0.847 ... ..$ sigma.theta: num [1:30000] 9.056 8.592 11.56 0.645 5.428 ... ..$ deviance : num [1:30000] 61 60.7 69.2 59.4 60.6 ... $ sims.matrix : num [1:30000, 1:11] 5.65 12.91 9.5 6.81 5.63 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:11] "theta[1]" "theta[2]" "theta[3]" "theta[4]" ... $ summary : num [1:11, 1:9] 11.75 8.08 6.32 7.75 5.47 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:11] "theta[1]" "theta[2]" "theta[3]" "theta[4]" ... .. ..$ : chr [1:9] "mean" "sd" "2.5%" "25%" ... $ mean :List of 4 ..$ theta : num [1:8(1d)] 11.75 8.08 6.32 7.75 5.47 ... ..$ mu.theta : num 8.14 ..$ sigma.theta: num 6.62 ..$ deviance : num 60.4 $ sd :List of 4 ..$ theta : num [1:8(1d)] 8.4 6.41 7.92 6.6 6.51 ... ..$ mu.theta : num 5.28 ..$ sigma.theta: num 5.68 ..$ deviance : num 2.21 $ median :List of 4 ..$ theta : num [1:8(1d)] 10.61 8.08 6.85 7.82 5.94 ... ..$ mu.theta : num 8.14 ..$ sigma.theta: num 5.25 ..$ deviance : num 60 $ root.short : chr [1:4] "theta" "mu.theta" "sigma.theta" "deviance" $ long.short :List of 4 ..$ : int [1:8] 4 5 6 7 8 9 10 11 ..$ : int 2 ..$ : int 3 ..$ : int 1 $ dimension.short: num [1:4] 1 0 0 0 $ indexes.short :List of 4 ..$ :List of 8 .. ..$ : num 1 .. ..$ : num 2 .. ..$ : num 3 .. ..$ : num 4 .. ..$ : num 5 .. ..$ : num 6 .. ..$ : num 7 .. ..$ : num 8 ..$ : NULL ..$ : NULL ..$ : NULL $ last.values :List of 3 ..$ :List of 3 .. ..$ theta : num [1:8] 25.58 10.95 11.31 24.44 4.85 ... .. ..$ mu.theta : num 8.74 .. ..$ sigma.theta: num 7.42 ..$ :List of 3 .. ..$ theta : num [1:8] 29.03 9.48 8.4 -1.63 2.23 ... .. ..$ mu.theta : num 6.8 .. ..$ sigma.theta: num 13.2 ..$ :List of 3 .. ..$ theta : num [1:8] 1.255 -3.145 0.201 0.803 -0.151 ... .. ..$ mu.theta : num 0.648 .. ..$ sigma.theta: num 1.49 $ isDIC : logi TRUE $ DICbyR : logi FALSE $ pD : num 63.2 $ DIC : num 2.8 $ model.file : chr "/Library/Frameworks/R.framework/Versions/2.13/Resources/library/R2WinBUGS/model/schools.txt" - attr(*, "class")= chr "bugs"
コメント 0