R/geweke.R

In coda: Output Analysis and Diagnostics for MCMC

"geweke.diag" <-  function (x, frac1 = 0.1, frac2 = 0.5) 
{
    if (frac1 < 0 || frac1 > 1) {
        stop("frac1 invalid")
    }
    if (frac2 < 0 || frac2 > 1) {
        stop("frac2 invalid")
    }
    if (frac1 + frac2 > 1) {
        stop("start and end sequences are overlapping")
    }
    if (is.mcmc.list(x)) {
        return(lapply(x, geweke.diag, frac1, frac2))
    }
    x <- as.mcmc(x)
    xstart <- c(start(x), floor(end(x) - frac2 * (end(x) - start(x))))
    xend <- c(ceiling(start(x) + frac1 * (end(x) - start(x))), end(x))
    y.variance <- y.mean <- vector("list", 2)
    for (i in 1:2) {
        y <- window(x, start = xstart[i], end = xend[i])
        y.mean[[i]] <- apply(as.matrix(y), 2, mean)
        y.variance[[i]] <- spectrum0.ar(y)$spec/niter(y)
    }
    z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])
    out <- list(z = z, frac = c(frac1, frac2))
    class(out) <- "geweke.diag"
    return(out)
}

"geweke.plot" <-
  function (x, frac1 = 0.1, frac2 = 0.5, nbins = 20, 
            pvalue = 0.05, auto.layout = TRUE, ask, ...) 
{
  if (missing(ask)) {
    ask <- if (is.R()) {
      dev.interactive()
    }
    else {
      interactive()
    }
  }
    
  x <- as.mcmc.list(x)
  oldpar <- NULL
  on.exit(par(oldpar))
  if (auto.layout) 
    oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 
                    Nparms = nvar(x)))
  ystart <- seq(from = start(x), to = (start(x) + end(x))/2, length = nbins)
  if (is.R())
    gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 
               dimnames = list(ystart, varnames(x), chanames(x)))
  else
    gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 
               dimnames = list(ystart, varnames(x), chanames(x)))
               
  for (n in 1:length(ystart)) {
    geweke.out <- geweke.diag(window(x, start = ystart[n]), 
                              frac1 = frac1, frac2 = frac2)
    for (k in 1:nchain(x)) gcd[n, , k] <- geweke.out[[k]]$z
  }
  climit <- qnorm(1 - pvalue/2)
  for (k in 1:nchain(x)) for (j in 1:nvar(x)) {
    ylimit <- max(c(climit, abs(gcd[, j, k])))
    plot(ystart, gcd[, j, k], type = "p", xlab = "First iteration in segment", 
         ylab = "Z-score", pch = 4, ylim = c(-ylimit, ylimit), 
         ...)
    abline(h = c(climit, -climit), lty = 2)
    if (nchain(x) > 1) {
      title(main = paste(varnames(x, allow.null = FALSE)[j], 
              " (", chanames(x, allow.null = FALSE)[k], ")", 
              sep = ""))
    }
    else {
      title(main = paste(varnames(x, allow.null = FALSE)[j], 
              sep = ""))
    }
    if (k==1 && j==1)
       oldpar <- c(oldpar, par(ask = ask))
  }
  invisible(list(start.iter = ystart, z = gcd))
}

"print.geweke.diag" <- function (x, digits = min(4, .Options$digits), ...) 
  ## Print method for output from geweke.diag
{
  cat("\nFraction in 1st window =", x$frac[1])
  cat("\nFraction in 2nd window =", x$frac[2], "\n\n")
  print.default(x$z, digits = digits, ...)
  cat("\n")
  invisible(x)
}