1"geweke.diag" <- function (x, frac1 = 0.1, frac2 = 0.5) 2{ 3 if (frac1 < 0 || frac1 > 1) { 4 stop("frac1 invalid") 5 } 6 if (frac2 < 0 || frac2 > 1) { 7 stop("frac2 invalid") 8 } 9 if (frac1 + frac2 > 1) { 10 stop("start and end sequences are overlapping") 11 } 12 if (is.mcmc.list(x)) { 13 return(lapply(x, geweke.diag, frac1, frac2)) 14 } 15 x <- as.mcmc(x) 16 xstart <- c(start(x), floor(end(x) - frac2 * (end(x) - start(x)))) 17 xend <- c(ceiling(start(x) + frac1 * (end(x) - start(x))), end(x)) 18 y.variance <- y.mean <- vector("list", 2) 19 for (i in 1:2) { 20 y <- window(x, start = xstart[i], end = xend[i]) 21 y.mean[[i]] <- apply(as.matrix(y), 2, mean) 22 y.variance[[i]] <- spectrum0.ar(y)$spec/niter(y) 23 } 24 z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]]) 25 out <- list(z = z, frac = c(frac1, frac2)) 26 class(out) <- "geweke.diag" 27 return(out) 28} 29 30"geweke.plot" <- 31 function (x, frac1 = 0.1, frac2 = 0.5, nbins = 20, 32 pvalue = 0.05, auto.layout = TRUE, ask, ...) 33{ 34 if (missing(ask)) { 35 ask <- if (is.R()) { 36 dev.interactive() 37 } 38 else { 39 interactive() 40 } 41 } 42 43 x <- as.mcmc.list(x) 44 oldpar <- NULL 45 on.exit(par(oldpar)) 46 if (auto.layout) 47 oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 48 Nparms = nvar(x))) 49 ystart <- seq(from = start(x), to = (start(x) + end(x))/2, length = nbins) 50 if (is.R()) 51 gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 52 dimnames = list(ystart, varnames(x), chanames(x))) 53 else 54 gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 55 dimnames = list(ystart, varnames(x), chanames(x))) 56 57 for (n in 1:length(ystart)) { 58 geweke.out <- geweke.diag(window(x, start = ystart[n]), 59 frac1 = frac1, frac2 = frac2) 60 for (k in 1:nchain(x)) gcd[n, , k] <- geweke.out[[k]]$z 61 } 62 climit <- qnorm(1 - pvalue/2) 63 for (k in 1:nchain(x)) for (j in 1:nvar(x)) { 64 ylimit <- max(c(climit, abs(gcd[, j, k]))) 65 plot(ystart, gcd[, j, k], type = "p", xlab = "First iteration in segment", 66 ylab = "Z-score", pch = 4, ylim = c(-ylimit, ylimit), 67 ...) 68 abline(h = c(climit, -climit), lty = 2) 69 if (nchain(x) > 1) { 70 title(main = paste(varnames(x, allow.null = FALSE)[j], 71 " (", chanames(x, allow.null = FALSE)[k], ")", 72 sep = "")) 73 } 74 else { 75 title(main = paste(varnames(x, allow.null = FALSE)[j], 76 sep = "")) 77 } 78 if (k==1 && j==1) 79 oldpar <- c(oldpar, par(ask = ask)) 80 } 81 invisible(list(start.iter = ystart, z = gcd)) 82} 83 84"print.geweke.diag" <- function (x, digits = min(4, .Options$digits), ...) 85 ## Print method for output from geweke.diag 86{ 87 cat("\nFraction in 1st window =", x$frac[1]) 88 cat("\nFraction in 2nd window =", x$frac[2], "\n\n") 89 print.default(x$z, digits = digits, ...) 90 cat("\n") 91 invisible(x) 92} 93 94 95 96 97 98 99 100 101 102 103 104