Nothing
print.bugs <- function(x, digits.summary = 1, ...)
{
fround <- function(x, digits) { format(round(x, digits), nsmall=digits) }
if(!is.null(x$model.file))
cat("Inference for Bugs model at \"", x$model.file, "\", ", sep="")
if(!is.null(x$program))
cat("fit using ", x$program, ",", sep="")
cat("\n ", x$n.chains, " chains, each with ", x$n.iter,
" iterations (first ", x$n.burnin, " discarded)", sep = "")
if(x$n.thin > 1) cat(", n.thin =", x$n.thin)
cat("\n n.sims =", x$n.sims, "iterations saved\n")
print(round(x$summary, digits.summary), ...)
if(x$n.chains > 1) {
cat("\nFor each parameter, n.eff is a crude measure of effective sample size,")
cat("\nand Rhat is the potential scale reduction factor (at convergence, Rhat=1).\n")
}
if(x$isDIC) {
msgDICRule <- ifelse(x$DICbyR,
"(using the rule, pD = var(deviance)/2)", ## Gelman tweak
"(using the rule, pD = Dbar-Dhat)") ## BUGS
cat(paste("\nDIC info ", msgDICRule, "\n", sep=""))
if(length(x$DIC) == 1) {
cat("pD =", fround(x$pD, 1), "and DIC =", fround(x$DIC, 1))
} else if(length(x$DIC)>1) {
print(round(x$DIC, 1))
}
cat("\nDIC is an estimate of expected predictive error (lower deviance is better).\n")
}
invisible(x)
}
plot.bugs <- function (x, display.parallel = FALSE, ...){
oldpar <- par(no.readonly = TRUE) # code line i
on.exit(par(oldpar)) # code line i + 1
mar.old <- par("mar")
pty.old <- par(pty = "m")
mfrow.old <- par("mfrow")
if (TRUE)
layout(matrix(c(1,2),1,2))
else
par(mfrow = c(1,2))
bugs.plot.summary (x, ...)
bugs.plot.inferences (x, display.parallel, ...)
header <- ""
if(!is.null(x$model.file))
header <- paste(header, "Bugs model at \"", x$model.file, "\", ", sep="")
if(!is.null(x$program))
header <- paste(header, "fit using ", x$program, ", ", sep="")
header <- paste(header, x$n.chains, " chains, each with ",
x$n.iter, " iterations (first ", x$n.burnin, " discarded)", sep = "")
mtext(header, outer = TRUE, line = -1, cex = 0.7)
if (TRUE) par(pty = pty.old[[1]], mar = mar.old, mfrow = mfrow.old)
else invisible(par(pty = pty.old[[1]], mar = mar.old, mfrow = mfrow.old))
}
"bugs.plot.summary" <-
function (sims, ...){
isDIC <- sims$isDIC
oldpar <- par(no.readonly = TRUE) # code line i
on.exit(par(oldpar)) # code line i + 1
if (.Device=="windows" ||
(.Device=="null device" && options("device")=="windows")){
cex.names <- .7
cex.top <- .7
cex.points <- .7
max.length <- 50
min.width <- .01
}
else {
cex.names <- .7
cex.top <- .7
cex.points <- .3
max.length <- 80
min.width <- .005
}
summ <- sims$summary
sims.array <- sims$sims.array
n.chains <- sims$n.chains
n.parameters <- nrow(summ)
J0 <- unlist(lapply(sims$long.short, length))
if (isDIC) J0 <- J0[1:(length(J0)-1)] # don't display deviance summaries
J <- J0
total <- ceiling(sum(J+.5))
while ((total > max.length) && max(J)>1){### vielleicht optimieren ...
J[J==max(J)] <- max(J)-1
total <- ceiling(sum(J+.5))
}
if (TRUE){
pos <- -1
} else {
pos <- -1.5
}
ypos <- NULL
id <- NULL
ystart <- NULL
jj <- 1:J[1]
n.roots <- length(sims$root.short)
if (isDIC) n.roots <- n.roots-1 # don't display deviance summaries
ystart <- numeric(n.roots)
for (k in 1:n.roots){
ystart[k] <- pos
if (TRUE) {
ypos <- c(ypos, pos - seq(0, J[k]-1))
} else {
# In S-PLUS, increase the vertical spacing
ypos <- c(ypos, pos - 1.5*seq(0, J[k]-1))
}
id <- c(id, 1:J[k])
if (TRUE) {
pos <- pos - J[k] -.5
} else {
pos <- pos - 1.5*J[k] -0.75
}
if (k>1) jj <- c(jj, sum(J0[1:(k-1)]) + (1:J[k]))
}
if (TRUE){
bottom <- min(ypos)-1
} else {
bottom <- min(ypos)-1.5
}
med <- numeric(sum(J))
i80 <- matrix( , sum(J), 2)
i80.chains <- array (NA, c(sum(J), n.chains, 2))
for (j in 1:sum(J)){
med[j] <- median (sims.array[,,jj[j]])
i80[j,] <- quantile (sims.array[,,jj[j]], c(.1,.9))
for (m in 1:n.chains)
i80.chains[j,m,] <- quantile (sims.array[,m,jj[j]], c(.1,.9))
}
rng <- range (i80, i80.chains)
p.rng <- pretty(rng, n = 2)
b <- 2 / (max(p.rng) - min(p.rng))
a <- -b * p.rng[1]
par (mar=c(0,0,1,3))
plot (c(0,1), c(min(bottom, -max.length)-3,2.5),
ann=FALSE, bty="n", xaxt="n", yaxt="n", type="n")
W <- max(strwidth(unlist(dimnames(summ)[[1]]), cex=cex.names))
B <- (1-W)/3.6
A <- 1-3.5*B
B <- (1-A)/3.5
b <- B*b
a <- A + B*a
text (A+B*1, 2.5, "80% interval for each chain", cex=cex.top)
lines (A+B*c(0,2), c(0,0))
lines (A+B*c(0,2), rep(bottom,2))
if(n.chains > 1){
text (A+B*3, 2.6, "R-hat", cex=cex.top)
lines (A+B*c(2.5,3.5), c(0,0))
lines (A+B*c(2.5,3.5), rep(bottom,2))
}
#
# line at zero
#
if (min(p.rng)<0 & max(p.rng)>0)
lines (rep(a,2), c(0,bottom), lwd=.5, col="gray")
for (x in p.rng){
text (a+b*x, 1, x, cex=cex.names)
lines (rep(a+b*x,2), c(0,-.2))
text (a+b*x, bottom-1, x, cex=cex.names)
lines (rep(a+b*x,2), bottom+c(0,.2))
}
if(n.chains > 1)
for (x in seq(1,2,.5)){
text (A+B*(1.5+seq(1,2,.5)), rep(1,3), c("1","1.5","2+"), cex=cex.names)
lines (A+B*rep(1.5+x,2), c(0,-.2))
text (A+B*(1.5+seq(1,2,.5)), rep(bottom-1,3), c("1","1.5","2+"),
cex=cex.names)
lines (A+B*rep(1.5+x,2), bottom+c(0,.2))
}
for (j in 1:sum(J)){
name <- dimnames(summ)[[1]][jj[j]]
if (id[j]==1)
if (TRUE) {
text (0, ypos[j], name, adj=0, cex=cex.names)
} else {
# in S-PLUS, strwidth is an upper bound on the length of the string,
# so we must align the brackets differently than in R
pos <- as.vector(regexpr("[[]", name))
text (0, ypos[j], substring(name, 1, pos-1), adj=0, cex=cex.names)
text (strwidth(substring(name,1,pos-1),cex=cex.names),
ypos[j], substring(name, pos, nchar(name)), adj=0, cex=cex.names)
}
else {
pos <- as.vector(regexpr("[[]", name))
text (strwidth(substring(name,1,pos-1),cex=cex.names),
ypos[j], substring(name, pos, nchar(name)), adj=0, cex=cex.names)
}
for (m in 1:n.chains){
interval <- a + b*i80.chains[j,m,]
if (interval[2]-interval[1] < min.width)
interval <- mean(interval) + c(-1,1)*min.width/2
lines (interval, rep(ypos[j]-.1*(m-(n.chains+1)/2),2), lwd=1, col=m+1)
if(n.chains > 1)
points (A+B*(1.5 + min(max(summ[jj[j],"Rhat"],1),2)), ypos[j], pch=20, cex=cex.points)
}
}
for (k in 1:n.roots){
if (J[k]<J0[k]) text (-.015, ystart[k], "*", cex=cex.names,
col="red")
}
if (sum(J!=J0)>0) text (0, bottom-3,
"* array truncated for lack of space", adj=0, cex=cex.names, col="red")
}
"bugs.plot.inferences" <-
function (sims, display.parallel, ...){
oldpar <- par(no.readonly = TRUE) # code line i
on.exit(par(oldpar)) # code line i + 1
if (.Device=="windows" ||
(.Device=="null device" && options("device")=="windows")){
cex.names <- .7
cex.axis <- .6
cex.tiny <- .4
cex.points <- .7
standard.width <- 30
max.width <- 40
min.width <- .02
}
else {
cex.names <- .7
cex.axis <- .6
cex.tiny <- .4
cex.points <- .3
standard.width <- 30
max.width <- 40
min.width <- .01
}
rootnames <- sims$root.short
n.roots <- length(rootnames)
sims.array <- sims$sims.array
n.chains <- sims$n.chains
dimension.short <- sims$dimension.short
indexes.short <- sims$indexes.short
long.short <- sims$long.short
height <- .6
par (mar=c(0,0,1,0))
plot (c(0,1), c(-n.roots-.5,-.4),
ann=FALSE, bty="n", xaxt="n", yaxt="n", type="n")
W <- max(strwidth(rootnames, cex=cex.names))
B <- (1-W)/3.8
A <- 1-3.5*B
if (display.parallel)
text (A, -.4, "80% interval for each chain", adj=0, cex=cex.names)
else
text (A, -.4, "medians and 80% intervals", adj=0, cex=cex.names)
num.height <- strheight (1:9, cex=cex.tiny)
for (k in 1:n.roots){
text (0, -k, rootnames[k], adj=0, cex=cex.names)
J <- min (length(long.short[[k]]), max.width)
if (k==1)
index <- 1:J
else
index <- sum (unlist(lapply(long.short,length))[1:(k-1)]) + 1:J
spacing <- 3.5/max(J,standard.width)
med <- numeric(J)
i80 <- matrix( , J, 2)
med.chains <- matrix( , J, sims$n.chains)
i80.chains <- array(NA, c(J, sims$n.chains, 2))
for (j in 1:J){
med[j] <- median(sims.array[,,index[j]])
i80[j,] <- quantile(sims.array[,,index[j]], c(.1,.9))
for (m in 1:n.chains){
med.chains[j,m] <- quantile (sims.array[,m,index[j]], .5)
i80.chains[j,m,] <- quantile (sims.array[,m,index[j]], c(.1,.9))
}
}
rng <- range (i80, i80.chains)
p.rng <- pretty(rng, n = 2)
b <- height/(max(p.rng) - min(p.rng))
a <- -(k+height/2) - b*p.rng[1]
lines (A+c(0,0), -k+c(-height/2,height/2))
#
# plot a line at zero (if zero is in the range of the mini-plot)
#
if (min(p.rng)<0 & max(p.rng)>0)
lines (A+B*spacing*c(0,J+1), rep (a,2), lwd=.5, col="gray")
for (x in p.rng){
text (A-B*.2, a+b*x, x, cex=cex.axis)
lines (A+B*c(-.05,0), rep(a+b*x,2))
}
for (j in 1:J){
if (display.parallel){
for (m in 1:n.chains){
interval <- a + b*i80.chains[j,m,]
if (interval[2]-interval[1] < min.width)
interval <- mean(interval) + c(-1,1)*min.width/2
lines (A+B*spacing*rep(j+.6*(m-(n.chains+1)/2)/n.chains,2),
interval, lwd=.5, col=m+1)
}
}
else {
lines (A+B*spacing*rep(j,2), a + b*i80[j,], lwd=.5)
for (m in 1:n.chains)
# points (A+B*spacing*j, a + b*med[j], pch=20, cex=cex.points)
points (A+B*spacing*j, a + b*med.chains[j,m], pch=20, cex=cex.points,
col=m+1)
}
dk <- dimension.short[k]
if (dk>0){
for (m in 1:dk){
index0 <- indexes.short[[k]][[j]][m]
if (j==1)
text(A+B*spacing*j, -k-height/2-.05-num.height*(m-1), index0,
cex=cex.tiny)
else if (index0!=indexes.short[[k]][[j-1]][m] &
(index0%%(floor(log10(index0)+1))==0))
text(A+B*spacing*j, -k-height/2-.05-num.height*(m-1), index0,
cex=cex.tiny)
}
}
}
if (J<length(long.short[[k]])) text (-.015, -k, "*",
cex=cex.names, col="red")
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.