"acf" <-
function (x, lag.max, type = "correlation", plot = T) 
{
        legal.types <- c("correlation", "covariance")
        check.type <- pmatch(type, legal.types)
        if (is.na(check.type)) 
                stop("Unknown type")
        else if (check.type == 0) 
                stop("Ambiguous type")
        else type <- legal.types[check.type]
        series <- deparse(substitute(x))
        x.freq <- frequency(as.ts(x))
        attr(x, "tsp") <- NULL
        x <- as.matrix(x)
        n.used <- nrow(x)
        if (missing(lag.max)) 
                lag.max <- floor(10 * log10(n.used))
        lag.max <- min(lag.max, n.used - 1)
        lag <- seq(from = 0, to = lag.max)
        acf <- matrix(nrow = lag.max + 1, ncol = ncol(x))
        x1 <- x2 <- sweep(x, 2, apply(x, 2, mean))
        for (i in 1:nrow(acf)) {
                acf[i, ] <- apply(x1 * x2, 2, sum)/n.used
                x1 <- x1[-1, , drop = F]
                x2 <- x2[-nrow(x2), , drop = F]
        }
        if (type == "correlation") 
                acf <- sweep(acf, 2, acf[1, ], "/")
        if (ncol(acf) == 1) 
                acf <- as.vector(acf)
        return(acf = acf, lag = lag/x.freq, n.used = n.used, 
                type = type, series = series)
}
"as.mcmc" <-
function (x) 
if (is.mcmc(x)) x else mcmc(x)
"as.ts.mcmc" <-
function (x) 
{
        x <- drop(x)
        if (length(dim(x)) == 3) 
                stop("Can't coerce to time series")
        else ts(x, start = start(x), end = end(x), deltat = thin(x))
}
"autocorr" <-
function (x, lags = c(1, 5, 10, 50) * thin(x)) 
{
        y <- array(dim = c(nvar(x), length(lags), nchain(x)))
        dimnames(y) <- list(varnames(x), paste("Lag", lags), 
                chanames(x))
        for (i in 1:nchain(x)) {
                acf.out <- acf(as.ts.mcmc(x[, , i]), lag.max = max(lags), 
                        plot = F)$acf
                y[, , i] <- if (is.matrix(acf.out)) 
                        t(acf.out[lags%/%thin(x) + 1, ])
                else acf.out[lags%/%thin(x) + 1]
        }
        y <- drop(y)
        return(y)
}
"autocorr.plot" <-
function (x, lag.max, auto.layout = T) 
{
        oldpar <- NULL
        if (auto.layout) 
                oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 
                        Nparms = nvar(x)))
        oldpar <- c(oldpar, set.scale())
        on.exit(par(oldpar))
        title.scale <- par("cex") * c(1, 0.85, 0.7)[coda.options("ps.plot")]
        scale <- par("cex") * c(1, 0.8, 0.6)[coda.options("ps.plot")]
        for (i in 1:nchain(x)) {
                xacf <- if (missing(lag.max)) 
                        acf(as.ts.mcmc(x[, , i]), plot = F)
                else acf(as.ts.mcmc(x[, , i]), lag.max = lag.max, 
                        plot = F)
                for (j in 1:nvar(x)) {
                        plot(xacf$lag, xacf$acf[, j], type = "h", 
                                ylab = "Autocorrelation", xlab = "Lag", 
                                cex = scale, ylim = c(-1, 1))
                        title(paste(varnames(x)[j], ":", chanames(x)[i], 
                                sep = ""), cex = title.scale)
                        if (j == nvar(x) && i == nchain(x)) 
                                stop.plot <- FALSE
                        else if (stop.plot <- mpause2()) 
                                break
                }
                par(new = F)
                if (stop.plot) 
                        break
        }
        invisible(x)
}
"change.tfoption" <-
function (string, option) 
{
        current.value <- coda.options(option)
        if (!is.logical(current.value)) 
                stop("Invalid option: must take logical values")
        wrd <- ifelse(current.value, " (Y/n)?\n:", " (y/N)?\n:")
        cat("\n", string, wrd, sep = "")
        ans <- readline()
        changeit <- (current.value && (ans == "N" | ans == "n")) || 
                (!current.value && (ans == "Y" | ans == "y"))
        if (changeit) {
                arg <- list(!current.value)
                names(arg) <- option
                coda.options(arg)
        }
        return()
}
"coda.credits" <-
function (file = "") 
{
        credits <- c("  _______________________________________________________________\n", 
                "|                                                               |\n", 
                "|                 Welcome to CODA for R!                        |\n", 
                "|  Convergence Diagnostics and Output Analysis for BUGS output  |\n", 
                "|_______________________________________________________________|\n", 
                "|                                                               |\n", 
                "|                                                               |\n", 
                "|  Authors : Martyn Plummer, Nicky Best, Kate Cowles,           |\n",
		"|            Karen Vines                                        |\n", 
                "|                                                               |\n", 
                "|  CODA version 0.3 Copyright (c) 1995-8 MRC Biostatistics Unit |\n", 
                "|  This is free software and commes with ABSOLUTELY NO WARRANTY |\n", 
                "|_______________________________________________________________|\n", 
                "\n")
        cat(credits, file = file)
}
"coda.options" <-
function (...) 
{
        #Set and display coda options
        #Works like options() and par(), ie
        #MTP start
        #Displays current values if no arguments are given
        #Displays current values of selected arguments if arguments are names
        #Resets values if argument is a named list
        #Like par it returns a value (rather than a list) if asked to
        #display one argument 
        #In addition
        #Resets to defaults if the "default=T" argument is given.
        #TODO put in value checking. Currently only checks mode
        #MTP finish
        single <- FALSE
        if (nargs() == 0) {
                return(.Coda.Options)
        }
        else {
                args <- list(...)
                if (length(args) == 1) {
                        if (is.list(args[[1]])) 
                                args <- args[[1]]
                        else if (is.null(names(args))) 
                                single <- TRUE
                }
        }
        if (is.null(names(args))) {
                #Display options
                args <- unlist(args)
                value <- vector("list", length(args))
                names(value) <- args
                for (v in args) if (any(v == names(.Coda.Options))) 
                        value[v] <- .Coda.Options[v]
                if (single) 
                        return(value[[1]])
                else return(value)
        }
        else {
                #Set options
                # MTP if there is no copy of .Coda.Options in the 
                # global environment we need to make one, otherwise
                # assignment to list members fails
                if (!exists(".Coda.Options", frame = 1)) 
                        .Coda.Options <<- .Coda.Options
                oldvalue <- vector("list", length(args))
                names(oldvalue) <- names(args)
                if (any(names(args) == "default") && args$default == 
                        T) 
                        .Coda.Options <<- .Coda.Options.Default
                for (v in names(args)) if (any(v == names(.Coda.Options))) {
                        oldvalue[v] <- .Coda.Options[v]
                        if (is.null(args[[v]])) 
                                .Coda.Options[v] <<- list(NULL)
                        else if (mode(.Coda.Options[[v]]) == 
                                mode(args[[v]])) 
                                .Coda.Options[v] <<- args[v]
                }
                invisible(oldvalue)
        }
}
.Coda.Options.Default <- list(
        trace = T,
        densplot = T,
        lowess = F,
        sepplot = F,
        onepage = F,
        bandwidth = function(x){1.06 * min(sd(x), IQR(x)/1.34) *
                                length(x[!is.na(x)])^-0.2},
        pos.parms = 0,
        pos.parmsmat = 0,
        prop.parms = 0,
        prop.parmsmat = 0,
        batch.size = 25,
        digits = 3,
        quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975),
        frac1 = 0.1,
        frac2 = 0.5,
        q = 0.025,
        r = 0.005,
        s = 0.95,
        combine.stats = F,
        combine.corr = F,
        ps.plot = 1,
        halfwidth = 0.1,
        user.layout = F,
        mrows = 1,
        mcols = 1,
        gr.bin = 10,
        geweke.bin = 10,
        gr.max = 50,
        geweke.max = 50)

.Coda.Options <- .Coda.Options.Default
"codamenu" <-
function () 
{
        #
        # codamenu -- main menu driver for CODA system
        #
        # Original author: Kate Cowles (University of Nebraska Medical Centre, USA)
        # Modified by: Nicky Best (MRC Biostatistics Unit, Cambridge, UK)
        # Additional contributions: Karen Vines (MRC Biostatistics Unit, Cambridge, UK)
        # R support: Martyn Plummer (IARC, France)
        # 
        # Date of last update: 9/3/98
        #
        options(object.size = 1e+09)
        on.exit(tidy.up())
        file.menu <- c("Begin a new CODA session using BUGS output files", 
                "Begin a new CODA session using data saved from a previous CODA session", 
                "Quit")
        coda.credits()
	cat("CODA startup menu\n\n")
        pick <- menu(file.menu)
        if (pick == 0 || pick == 3) 
		return(invisible())
	else if (pick == 1)
		coda.dat <<- read.bugs.interactive()
	else if (pick == 2) {
                msg <- "\nEnter name of R data object saved from a previous CODA session:\n"
                while (T) {
                        cat(msg, "\n")
                        outname <- scan(what = character(), sep = "\n", 
                                nmax = 1, strip.white = T)
                        if (length(outname) == 0) 
                                msg <- "You must enter something"
                        else if (outname == "0") {
                                return()
                        }
                        else if (!exists(outname)) 
                                msg <- "Can\'t find this object"
                        else if (!is.mcmc(eval(parse(text = outname)))) 
                                msg <- "Not an mcmc object"
                        else {
                                coda.dat <<- eval(parse(text = outname))
                                break
                        }
                        msg <- paste(msg, "Please re-enter another name or type 0 to quit:\n", 
                                sep = "\n")
		}
        }
	else stop("Invalid option")
	#
	# Create working data (Isn't this a bit wasteful? This could be the third copy of the
	# data in the workspace... MTP)
	#
	work.dat <<- coda.dat
	#
        # Enter menu loop.
        #
        current.menu <- "codamenu.main"
        repeat {
                next.menu <- do.call(current.menu, NULL)
                if (next.menu == "quit") {
                        cat("Are you sure you want to quit? (y/N)\n")
                        if (any(readline() == c("y", "Y"))) 
                                break
                }
                else current.menu <- next.menu
        }
        invisible()
}
"codamenu.anal" <-
function () 
{
        # 
        # outanal -- Output analysis menu for CODA system 
        # 
        # Author: Kate Cowles 
        # 
        next.menu <- "codamenu.anal"
        mtitle <- "CODA Output Analysis Menu"
        cat("", mtitle, paste(rep("*", nchar(mtitle)), collapse = ""), 
                sep = "\n")
        choices <- c("Plots", "Statistics", "List/Change Options", 
                "Return to Main Menu")
        next.menu.list <- c("plots", "summary", "codamenu.options", 
                "codamenu.main")
        pick <- menu(choices)
        if (pick == 0) 
                next.menu <- "quit"
        else if (next.menu.list[pick] == "summary") 
                print(summary(work.dat, quantiles = coda.options("quantiles"), 
                        combine.chains = coda.options("combine.stats")), 
                        digits = coda.options("digits"))
        else if (next.menu.list[pick] == "plots") {
                plot(work.dat, trace = coda.options("trace"), 
                        density = coda.options("densplot"), smooth = coda.options("lowess"), 
                        auto.layout = !coda.options("user.layout"), 
                        bwf = coda.options("bandwidth"), combine.chains = !coda.options("sepplot"), 
                        one.page = coda.options("onepage"))
                mpause3()
        }
        else next.menu <- next.menu.list[pick]
        return(next.menu)
}
"codamenu.diags" <-
function () 
{
        #
        # diags -- Diagnostics menu for CODA system
        #
        # Author: Kate Cowles
        #
        next.menu <- "diags"
        while (next.menu == "diags") {
                cat("\nCODA Diagnostics Menu\n*********************\n")
                choices <- c("Geweke", "Gelman and Rubin", "Raftery and Lewis", 
                        "Heidelberger and Welch", "Autocorrelations", 
                        "Cross-Correlations", "List/Change Options", 
                        "Return to Main Menu")
                next.menu.list <- c("codamenu.diags.geweke", 
                        "codamenu.diags.gelman", "codamenu.diags.raftery", 
                        "codamenu.diags.heidel", "codamenu.diags.autocorr", 
                        "codamenu.diags.crosscorr", "codamenu.options", 
                        "codamenu.main")
                pick <- menu(choices)
                if (pick == 0) 
                        return("quit")
                else next.menu <- next.menu.list[pick]
        }
        return(next.menu)
}
"codamenu.diags.autocorr" <-
function () 
{
        #
        # autocorr -- Calculate Autocorrelations for CODA system
        #
        # Author: Kate Cowles
        # Modified by: Nicky Best
        #
        next.menu <- "codamenu.diags"
        cat("\nLAGS AND AUTOCORRELATIONS WITHIN EACH CHAIN:")
        cat("\n============================================\n\n")
        print(autocorr(work.dat), digits = coda.options("digits"))
        cat("\nAutocorrelation Plots Menu\n")
        cat("**************************\n")
        choices <- c("Plot autocorrelations", "Return to Diagnostics Menu")
        pick <- menu(choices)
        if (pick == 0) 
                next.menu <- "quit"
        else if (pick == 1) {
                autocorr.plot(work.dat, auto.layout = !coda.options("user.layout"))
                mpause3()
        }
        return(next.menu)
}
"codamenu.diags.crosscorr" <-
function () 
{
        next.menu <- "codamenu.diags.crosscorr"
        crosscorr.out <- crosscorr(work.dat, combine.chains = coda.options("combine.corr"))
        if (coda.options("combine.corr") & nchain(work.dat) > 
                1) 
                cat("Pooling over chains:", chanames(work.dat), 
                        sep = "\n", collapse = "\n")
        print(crosscorr.out, digits = coda.options("digits"))
        cat("\nCross Correlation Plots Menu\n")
        cat("****************************\n")
        choices <- c("Change options", "Plot Cross Correlations", 
                "Return to Diagnostics Menu")
        pick <- menu(choices)
        if (pick == 0) 
                next.menu <- "quit"
        else switch(pick, change.tfoption("Combine chains", "combine.corr"), 
                {
                        crosscorr.plot(work.dat, combine.chains = coda.options("combine.corr"))
                        mpause3()
                }, next.menu <- "codamenu.diags")
        return(next.menu)
}
"codamenu.diags.gelman" <-
function () 
{
        next.menu <- this.menu <- "codamenu.diags.gelman"
        if (nchain(work.dat) == 1) {
                cat("\nError: Gelman and Rubin convergence diagnostic requires more than one chain.\n\n")
                return("codamenu.diags")
        }
        z <- window(work.dat, start = niter(work.dat)/2)
        tol <- 1e-08
        if (any(apply(z[, , 1] - z[, , -1], c(2, 3), var) < tol)) {
                cat("\nError: 2nd half of two chains are identical for at least one variable\nGelman and Rubin convergence diagnostic requires 2 or more different chains.\n\n")
                next.menu <- "codamenu.diags"
        }
        else if (niter(work.dat) <= 50) {
                cat("\nError: you must have more than 50 iterations in the working data\n")
                next.menu <- "codamenu.diags"
        }
        else {
                codamenu.output.header("GELMAN AND RUBIN DIAGNOSTIC")
                print(gelman.diag(work.dat, transform = T), digits = coda.options("digits"))
        }
        while (next.menu == "codamenu.diags.gelman") {
                cat("\nGelman & Rubin Menu\n")
                cat("*************************\n")
                choices <- c("Trace Plots", "Shrink Factor Plots", 
                        "Change bin size for shrink plot", "Return to Diagnostics Menu")
                action.list <- c("TracePlot", "ShrinkPlot", "ChangeBin", 
                        "Return")
                pick <- menu(choices)
                if (pick == 0) 
                        next.menu <- "quit"
                else switch(action.list[pick], TracePlot = {
                        gelman.trace.plot(work.dat, auto.layout = !coda.options("user.layout"))
                        mpause3()
                }, ShrinkPlot = {
                        if (is.numeric(coda.options("gr.bin")) & 
                                niter(work.dat) <= 0 + coda.options("gr.bin")) {
                                cat("\n******* Error: *******")
                                cat("\nNot enough iterations available for plot\n")
                        }
                        else {
                                if (is.numeric(coda.options("gr.max")) & 
                                 niter(work.dat) < 50 + coda.options("gr.max")) {
                                 cat("\n******* Error: *******")
                                 cat("\nNot enough iterations available for plot\n")
                                }
                                else {
                                 max.bins <- coda.options("gr.max")
                                 bin.width <- coda.options("gr.bin")
                                 auto.layout <- !coda.options("user.layout")
                                 if (max.bins != 0 & bin.width != 
                                  0) 
                                  gelman.shrink.plot(work.dat, 
                                   max.bins = max.bins, bin.width = bin.width, 
                                   auto.layout = auto.layout)
                                 else if (max.bins != 0) 
                                  gelman.shrink.plot(work.dat, 
                                   max.bins = max.bins, auto.layout = auto.layout)
                                 else if (bin.width != 0) 
                                  gelman.shrink.plot(work.dat, 
                                   bin.width = bin.width, auto.layout = auto.layout)
                                 else gelman.shrink.plot(work.dat, 
                                  auto.layout = auto.layout)
                                 mpause3()
                                }
                        }
                }, ChangeBin = {
                        next.menu <- codamenu.options.gelman(this.menu)
                }, Return = {
                        next.menu <- "codamenu.diags"
                })
        }
        return(next.menu)
}
"codamenu.diags.geweke" <-
function () 
{
        this.menu <- next.menu <- "codamenu.diags.geweke"
        codamenu.output.header("GEWEKE CONVERGENCE DIAGNOSTIC (Z-score)")
        print(geweke.diag(work.dat, frac1 = coda.options("frac1"), 
                frac2 = coda.options("frac2")), digits = coda.options("digits"))
        cat("\nGeweke Plots Menu\n")
        cat("*****************\n")
        choices <- c("Change window size", "Plot Z-scores", "Change bin size for plot", 
                "Return to Diagnostics Menu")
        action.list <- c("ChangeWindow", "Plot", "ChangeBin", 
                "Return")
        while (next.menu == "codamenu.diags.geweke") {
                pick <- menu(choices)
                if (pick == 0) 
                        return("quit")
                switch(action.list[pick], ChangeWindow = {
                        next.menu <- codamenu.options.geweke.win(this.menu)
                        codamenu.output.header("")
                        print(geweke.diag(work.dat, frac1 = coda.options("frac1"), 
                                frac2 = coda.options("frac2")), 
                                digits = coda.options("digits"))
                }, Plot = {
                        if ((is.numeric(coda.options("geweke.bin")) && 
                                nrow(work.dat) < 50 + coda.options("geweke.bin")) || 
                                (is.numeric(coda.options("geweke.max")) && 
                                 nrow(work.dat) < 50 + coda.options("geweke.max"))) {
                                cat("\n******* Error: *******")
                                cat("\nNot enough iterations available for plot\n")
                        }
                        else {
                                frac1 <- coda.options("frac1")
                                frac2 <- coda.options("frac2")
                                max.bins <- coda.options("geweke.max")
                                bin.width <- coda.options("geweke.bin")
                                auto <- !coda.options("user.layout")
                                if (max.bins != 0 & bin.width != 
                                 0) 
                                 geweke.plot(work.dat, frac1 = frac1, 
                                  frac2 = frac2, max.bins = max.bins, 
                                  bin.width = bin.width, auto.layout = auto)
                                else if (max.bins != 0) 
                                 geweke.plot(work.dat, frac1 = frac1, 
                                  frac2 = frac2, max.bins = max.bins, 
                                  auto.layout = auto)
                                else if (bin.width != 0) 
                                 geweke.plot(work.dat, frac1 = frac1, 
                                  frac2 = frac2, bin.width = bin.width, 
                                  auto.layout = auto)
                                else geweke.plot(work.dat, frac1 = frac1, 
                                 frac2 = frac2, auto.layout = auto)
                                mpause3()
                        }
                }, ChangeBin = {
                        next.menu <- codamenu.options.geweke.bin(this.menu)
                }, Return = {
                        next.menu <- "codamenu.diags"
                })
        }
        return(next.menu)
}
codamenu.diags.heidel <- function() {
	this.menu <- "codamenu.diags.heidel"
	next.menu <- "codamenu.diags"
        title <- "HEIDELBERGER AND WELCH STATIONARITY AND INTERVAL HALFWIDTH TESTS:"
	codamenu.output.header(title)
        cat("Precision of halfwidth test =", coda.options("halfwidth"), "\n\n")
	print(heidel.diag(work.dat, eps=coda.options("halfwidth")), 
		digits=coda.options("digits"))
	choices <- c("Change precision", "Return to diagnostics menu")
	pick <- menu(choices)
	if(pick==0)
		next.menu <- "quit"
	else if (pick==1)
		next.menu <- codamenu.options.heidel(this.menu)
	return(next.menu)
}
"codamenu.diags.raftery" <-
function () 
{
	this.menu <- "codamenu.diags.raftery"
        print("RAFTERY AND LEWIS CONVERGENCE DIAGNOSTIC:\n")
	print(raftery.diag(work.dat, q = coda.options("q"), r=coda.options("r"), s=coda.options("s")), digits=coda.options("digits"))
	choices <- c("Change parameters", "Return to diagnostics menu")
	pick <- menu(choices)
	next.menu <- if(pick==0)
		"quit"
	else if (pick==1)
		codamenu.options.raftery(this.menu)
	else
		"codamenu.diags"
	return(next.menu)
}
"codamenu.main" <-
function () 
{
        # 
        # codamenu.main -- CODA main menu 
        # 
        # Author: Kate Cowles 
        # Modified by: Nicky Best 
        # and by: Martyn Plummer 
        # 
        cat("\nCODA Main Menu\n**************\n")
        choices <- c("Output Analysis", "Diagnostics", "List/Change Options", 
                "Quit")
        next.menu.list <- c("codamenu.anal", "codamenu.diags", 
                "codamenu.options", "quit")
        pick <- menu(choices)
        if (pick == 0) 
                next.menu <- "quit"
        else next.menu <- next.menu.list[pick]
        return(next.menu)
}
"codamenu.options" <-
function () 
{
        #
        # codamenu.options-- Main options menu for CODA system
        #
        # Author: Nicky Best
        #
        next.menu <- "codamenu.options"
        choices <- c("List current options", "Data Options", 
                "Plot Options", "Summary Statistics Options", 
                "Diagnostics Options", "Output Analysis", "Diagnostics", 
                "Main Menu")
        action.list <- c("ListOptions", "codamenu.options.data", 
                "codamenu.options.plot", "codamenu.options.stats", 
                "codamenu.options.diag", "codamenu.anal", "codamenu.diags", 
                "codamenu.main")
        cat("\nCODA Main Options Menu\n***********************\n")
        pick <- menu(choices)
        if (pick == 0) 
                return("quit")
        if (action.list[pick] == "ListOptions") {
                print.coda.options(data = T, stats = T, plots = T, 
                        diags = T)
                next.menu <- "codamenu.options"
        }
        else next.menu <- action.list[pick]
        return(next.menu)
}
"codamenu.options.data" <-
function () 
{
        next.menu <- "codamenu.options.data"
        #
        work.vars <- varnames(work.dat)
        work.chains <- chanames(work.dat)
        work.start <- start(work.dat)
        work.end <- end(work.dat)
        work.thin <- thin(work.dat)
        #
        choices <- c("List current data options", "Select variables for analysis", 
                "Select chains for analysis", "Select iterations for analysis", 
                "Select thinning interval", "Return to main options menu")
        action.list <- c("ListDataOptions", "SelectVars", "SelectChains", 
                "SelectIters", "SelectThinInterval", "MainOptionsMenu")
        cat("\nCODA Data Options Menu\n***********************\n")
        pick <- menu(choices)
        if (pick == 0) 
                return("quit")
        switch(action.list[pick], ListDataOptions = {
                print.coda.options(data = T)
        }, SelectVars = {
                work.vars <- multi.menu(colnames(coda.dat), "Select variables for analysis", 
                        c("VARIABLE NUMBER", "VARIABLE NAME"), 
                        allow.zero = F)
        }, SelectChains = {
                work.chains <- multi.menu(dimnames(coda.dat)[[3]], 
                        "Select chains for analysis:", c("CHAIN NUMBER", 
                                "CHAIN NAME"), allow.zero = F)
        }, SelectIters = {
                cat("\nIterations available = ", start(coda.dat), 
                        ":", end(coda.dat), "\n", sep = "")
                work.start <- read.and.check("Enter iteration you wish to start at", 
                        lower = start(coda.dat), upper = end(coda.dat), 
                        default = start(work.dat))
                work.end <- read.and.check("Enter iteration you wish to end at", 
                        lower = work.start, upper = end(coda.dat), 
                        default = end(work.dat))
        }, SelectThinInterval = {
                cat("\nThinning interval of full data = ", thin(coda.dat), 
                        "\n", sep = "")
                work.thin <- read.and.check("Enter thinning interval:", 
                        lower = thin(coda.dat), default = thin(work.dat))
        }, MainOptionsMenu = {
                next.menu <- "codamenu.options"
        })
        if (action.list[pick] != "ListDataOptions" && action.list[pick] != 
                "MainOptionsMenu") {
                cat("Recreating working data...\n")
                work.dat <<- window(coda.dat[, work.vars, work.chains], 
                        start = work.start, end = work.end, thin = work.thin)
        }
        return(next.menu)
}
"codamenu.options.diag" <-
function () 
{
        next.menu <- this.menu <- "codamenu.options.diag"
        mtitle <- "CODA Diagnostics Options Menu"
        cat("", mtitle, paste(rep("*", nchar(mtitle)), collapse = ""), 
                sep = "\n")
        choices <- c("Display current diagnostic options", "Window sizes for Geweke's diagnostic", 
                "Bin size for plotting Geweke's diagnostic", 
                "Bin size for plotting Gelman & Rubin's diagnostic", 
                "Parameters for Raftery & Lewis' diagnostic", 
                "Halfwidth precision for Heidelberger & Welch's diagnostic", 
                "Combine chains to calculate correlation matrix", 
                "Return to main options menu")
        pick <- menu(choices)
        if (pick == 0) 
                return("quit")
        switch(pick, print.coda.options(diags = T), next.menu <- codamenu.options.geweke.win(this.menu), 
                next.menu <- codamenu.options.geweke.bin(this.menu), 
                next.menu <- codamenu.options.gelman(this.menu), 
                next.menu <- codamenu.options.raftery(this.menu), 
                next.menu <- codamenu.options.heidel(this.menu), 
                {
                        change.tfoption("Do you want to combine all chains to calculate correlation matrix", 
                                "combine.corr")
                }, next.menu <- "codamenu.options")
        return(next.menu)
}
"codamenu.options.gelman" <-
function (next.menu) 
{
        cat("\nOptions for defining bin size to plot Gelman-Rubin-Brooks diagnostic:\n")
        choices <- c("Default: bin width = 10; maximum number of bins = 50", 
                "User-specified bin width", "User-specified total number of bins")
        pick <- menu(choices)
        if (pick == 0) 
                return("quit")
        switch(pick, {
                coda.options(gr.max = 50)
                coda.options(gr.bin = 10)
        }, {
                coda.options(gr.max = 0)
                default <- if (coda.options("gr.bin") == 0) 
                        10
                else coda.options("gr.bin")
                msg <- "Enter required bin width:"
                coda.options(gr.bin = read.and.check(msg, lower = 1, 
                        upper = niter(work.dat) - 50, default = default))
        }, {
                coda.options(gr.bin = 0)
                default <- if (coda.options("gr.max") == 0) 
                        50
                else coda.options("gr.max")
                msg <- "Enter total number of bins required:"
                coda.options(gr.max = read.and.check(msg, lower = 1, 
                        upper = niter(work.dat) - 50, default = default))
        })
        return(next.menu)
}
"codamenu.options.geweke.bin" <-
function (next.menu) 
{
        cat("\nOptions for defining bin size to plot Geweke's diagnostic:\n")
        choices <- c("Default: bin width = 10; maximum number of bins = 50", 
                "User-specified bin width", "User-specified total number of bins")
        pick <- menu(choices)
        if (pick == 0) 
                return("quit")
        switch(pick, {
                coda.options(geweke.max = 50)
                coda.options(geweke.bin = 10)
        }, {
                coda.options(geweke.max = 0)
                default <- if (coda.options("geweke.bin") == 
                        0) 
                        10
                else coda.options("geweke.bin")
                msg <- "Enter required bin width:"
                coda.options(geweke.bin = read.and.check(msg, 
                        lower = 1, upper = niter(work.dat) - 
                                50, default = default))
        }, {
                coda.options(geweke.bin = 0)
                default <- if (coda.options("geweke.max") == 
                        0) 
                        50
                else coda.options("geweke.max")
                msg <- "Enter total number of bins required:"
                coda.options(geweke.max = read.and.check(msg, 
                        lower = 1, upper = niter(work.dat) - 
                                50, default = default))
        })
        return(next.menu)
}
"codamenu.options.geweke.win" <-
function (next.menu) 
{
        msg1 <- "Enter fraction of chain to include in 1st window of \nGeweke's diagnostic:"
        msg2 <- "Enter fraction of chain to include in 2nd window of \nGeweke's diagnostic:"
        ans1 <- ans2 <- 1
        while (ans1 + ans2 >= 1) {
                ans1 <- read.and.check(msg1, lower = 0, upper = 1, 
                        default = coda.options("frac1"))
                ans2 <- read.and.check(msg2, lower = 0, upper = 1, 
                        default = coda.options("frac2"))
                # Check that sum of fractions doesn't exceed 1.0
                msg1 <- if (ans1 + ans2 == 1) 
                        "Error: Sum of fractions in 1st and 2nd windows equals 1.0\nYou must leave a gap in the chain between the 2 windows\nPlease re-enter fraction of chain to include in 1st window:"
                else "Error: Sum of fractions in 1st and 2nd windows exceeds 1.0\nPlease re-enter fraction of chain to include in 1st window:"
                msg2 <- "Now re-enter fraction of chain to include in 2nd window:"
        }
        coda.options(frac1 = ans1, frac2 = ans2)
        return(next.menu)
}
"codamenu.options.heidel" <-
function (next.menu) 
{
        coda.options(halfwidth = read.and.check("Enter precision for halfwidth test", 
                default = coda.options("halfwidth")))
        return(next.menu)
}
"codamenu.options.plot" <-
function () 
{
        next.menu <- "codamenu.options.plot"
        choices <- c("Show current plotting options", "Plot trace of samples", 
                "Plot kernel density estimate", "Add smooth line through trace plot", 
                "Separate plots per chain", "Separate page of plots per variable", 
                "Specify page layout for plots", "Select bandwidth function for kernel smoothing", 
                "Options for saving plots to postscript file", 
                "Return to main options menu")
        cat("\nCODA Plotting Options Menu\n***************************\n")
        #
        pick <- menu(choices)
        # Scientific notation  
        #
        if (pick == 0) 
                return("quit")
        switch(pick, print.coda.options(plots = T), change.tfoption(choices[2], 
                "trace"), change.tfoption(choices[3], "densplot"), 
                change.tfoption(choices[4], "lowess"), change.tfoption(choices[5], 
                        "sepplot"), change.tfoption(choices[6], 
                        "onepage"), {
                        change.tfoption("Do you want to specify your own page layout for the plots", 
                                "user.layout")
                        if (coda.options("user.layout")) {
                                mrows <- read.and.check("Enter number of rows of plots per page (maximum=7)", 
                                 lower = 1, upper = 7)
                                mcols <- read.and.check("Enter number of columns of plots per page (maximum=8)", 
                                 lower = 1, upper = 8)
                                coda.options(mrows = mrows)
                                coda.options(mcols = mcols)
                                par(mfrow = c(mrows, mcols))
                        }
                }, {
                        next.menu <- "codamenu.options.plot.kernel"
                }, {
                        next.menu <- "codamenu.options.plot.ps"
                }, NULL)
        if (pick == length(choices)) 
                next.menu <- "codamenu.options"
        return(next.menu)
}
"codamenu.options.plot.kernel" <-
function () 
{
        if (!coda.options("densplot")) {
                cat("\nNo density plots requested - this option is irrelevant\n")
        }
        else {
                cat("\nSelect kernel bandwidth function required:\n")
                kernel.menu <- c("Smooth (0.25 * sample range)", 
                        "Coarse (Silverman 1986 eqn. 3.28 & 3.30)", 
                        "User-defined function", "Return to Plotting Options Menu")
                pick1 <- menu(kernel.menu)
                if (pick1 == 0) 
                        return("quit")
                switch(pick1, {
                        bw <- function(y) {
                                (max(y) - min(y))/4
                        }
                        coda.options(bandwidth = bw)
                }, {
                        bw <- function(x) {
                                1.06 * min(sd(x), IQR(x)/1.34) * 
                                 length(x)^-0.2
                        }
                        coda.options(bandwidth = bw)
                }, {
                        #
                        # Allow user to enter their own bandwidth function.
                        # This is read in as a series of character strings
                        # which are then collapsed into a single string with 
                        # no white space, before parsing and evaluating
                        # 
                        cat("DANGER! DANGER WILL ROBINSON!\n")
                        cat("If you make a mistake here you will crash\n")
                        cat("the codamenu session. Sorry, but there is\n")
                        cat("no way round this at the moment (R 0.61)\n")
                        cat("Do you want to try? (y/N)\n")
                        ans <- readline()
                        if (ans != "y" && ans != "Y") 
                                return()
                        func.OK <- F
                        while (!func.OK) {
                                cat("Enter bandwidth as a function of y, the sampled values, e.g. \n(max(y) - min(y)) / 4\n")
                                ans <- scan(what = character())
                                if (length(ans) > 0) {
                                 bw <- "function(y){"
                                 for (i in 1:length(ans)) {
                                  bw <- paste(bw, ans[i], sep = "")
                                 }
                                 bw <- paste(bw, "}", sep = "")
                                 #
                                 bw <- eval(parse(text = bw))
                                 # Carry out simple test to check whether the
                                 # function entered makes sense
                                 #
                                 func.OK <- test.bandwidth()
                                }
                        }
                        coda.options(bandwidth = bw)
                }, NULL)
        }
        return("codamenu.options.plot")
}
"codamenu.options.plot.ps" <-
function () 
{
        title <- "Select options required when saving plots to postscript files"
        cat("\n", title, "\n", sep = "")
        choices <- c("Portrait", "Landscape")
        pick <- menu(choices)
        if (pick == 0) 
                return("quit")
        else coda.options(ps.orientation = c("portrait", "landscape")[pick])
        if (.Device == "X11") 
                x11(orientation = coda.options("ps.orientation"))
        else if (.Device == "Win32") 
                windows(orientation = coda.options("ps.orientation"))
        return("codamenu.options.plot")
}
"codamenu.options.raftery" <-
function (next.menu) 
{
        coda.options(q = read.and.check("Enter quantile to be estimated:", 
                lower = 0, upper = 1, default = coda.options("q")))
        coda.options(r = read.and.check("Enter required precision:", 
                upper = coda.options("q"), default = coda.options("r")))
        coda.options(s = read.and.check("Enter required probability:", 
                lower = 0, upper = 1, default = coda.options("s")))
        return(next.menu)
}
"codamenu.options.stats" <-
function () 
{
        #
        # codamenu.options.stats -- Summary statistics options menu for CODA system
        #
        # Author: Nicky Best
        #
        #
        # Check that quantiles are between 0 and 1
        #
        next.menu <- "codamenu.options.stats"
        cat("\nCODA Summary Statistics Defaults Menu \n*************************************\n")
        choices <- c("Display current statistics options", "Combine chains for summary statistics", 
                "Quantiles for summary statistics", 
                "Number of significant digits for printing", 
                "Return to main options menu")
        pick <- menu(choices)
        if (pick == 0) 
                return("quit")
        switch(pick, print.coda.options(stats = T), {
                mssg <- "Do you want to combine all chains when calculating summary statistics"
                change.tfoption(mssg, "combine.stats")
        }, {
                mssg <- paste("Enter quantiles required, separated by commas\n(Default =", paste(coda.options("quantiles"), collapse=", "))
                repeat {
                        cat("\n", mssg, "\n")
                        ans <- as.numeric(scan(what = character(), 
                                sep = ",", quiet = T, nlines = 1))
                        if (any(is.na(ans))) 
                                mssg <- "You must enter numeric values"
                        else if (any(ans >= 1) || any(ans <= 
                                0)) 
                                mssg <- "You must enter values between 0 and 1"
                        else break
                }
                if (length(ans) > 0) 
                        coda.options(quantiles = sort(ans))
        }, {
                mssg <- "Enter number of significant digits to be printed"
                ans <- read.and.check(mssg, what = integer(), 
                        lower = 0, default=coda.options("digits"))
                coda.options(digits = ans)
        }, {
                next.menu <- "codamenu.options"
        })
        return(next.menu)
}
"codamenu.output.header" <-
function (title) 
{
        #
        # A short header: common to most codamenu output
        #
        cat("\n", title, sep = "")
        cat("\n", paste(rep("=", nchar(title)), collapse = ""), 
                "\n\n", sep = "")
        cat("Iterations used = ", start(work.dat), ":", end(work.dat), 
                "\n", sep = "")
        cat("Thinning interval =", thin(work.dat), "\n")
        cat("Sample size per chain =", niter(work.dat), "\n\n")
        return()
}
"crosscorr" <-
function (x, combine.chains = F) 
{
        x <- as.mcmc(x)
        Nvars <- dim(x)[2]
        Nchains <- dim(x)[3]
        if (combine.chains & Nchains > 1) 
                y <- cor(matrix(aperm(x, c(2, 1, 3)), ncol = nvar(x), 
                        byrow = T))
        else {
                y <- array(apply(x, 3, cor), dim = c(Nvars, Nvars, 
                        Nchains))
                dimnames(y)[[3]] <- dimnames(x)[[3]]
        }
        dimnames(y)[1:2] <- list(dimnames(x)[[2]], dimnames(x)[[2]])
        y <- drop(y)
        return(y)
}
"crosscorr.plot" <-
function (x, combine.chains = F, col = topo.colors(10), auto.layout = F) 
{
        x <- as.mcmc(x)
        Nvar <- nvar(x)
        pcorr <- crosscorr(x, combine.chains = combine.chains)
        if (is.matrix(pcorr)) 
                pcorr <- array(pcorr, dim = c(dim(pcorr), 1))
        dens <- ((pcorr + 1) * length(col))%/%2 + (pcorr < 1) + 
                (pcorr < -1)
        cutoffs <- format(seq(from = 1, to = -1, length = length(col) + 
                1), digits = 2)
        leg <- paste("(", cutoffs[-1], ",", cutoffs[-length(cutoffs)], 
                "]", sep = "")
        oldpar <- NULL
        if (auto.layout) 
                oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 
                        Nparms = nvar(x)))
        oldpar <- c(oldpar, set.scale())
        on.exit(par(oldpar))
        par(pty = "s")
        for (i in 1:dim(pcorr)[3]) {
                plot(0, 0, type = "n", xlim = c(0, Nvar), ylim = c(0, 
                        Nvar), xlab = "", ylab = "", xaxt = "n", 
                        yaxt = "n")
                axis(1, at = c(1:Nvar) - 0.5, labels = abbreviate(colnames(x), 
                        minlength = 7))
                axis(2, at = 1:Nvar - 0.5, labels = abbreviate(colnames(x), 
                        minlength = 7)[Nvar:1])
                for (cl in 1:Nvar) {
                        for (rw in 1:(Nvar - cl + 1)) polygon(y = c(cl - 
                                1, cl - 1, cl, cl, cl - 1), x = c(rw - 
                                1, rw, rw, rw - 1, rw - 1), col = col[dens[nrow(dens) - 
                                cl + 1, rw, i]])
                }
                oldpar <- c(oldpar, par(adj = 0.5))
                legend(Nvar, Nvar, legend = leg, fill = col[length(col):1], 
                        xjust = 1, cex = par("cex"), y.intersp = 1.2 * 
                                par("cex"))
                if (!combine.chains) 
                        title(dimnames(x)[[3]][i])
                if (i != dim(pcorr)[3] && mpause2()) 
                        break
        }
        if (combine.chains) 
                title("Combined")
}
"densplot" <-
function (x, combine.chains = F, bwf) 
{
        x <- as.mcmc(x)
        do.densplot <- function(j, i) {
                #NB Using R scoping rules
                #
                if (missing(i)) 
                        y <- as.vector(x[, j, , drop = T])
                else y <- as.vector(x[, j, i, drop = T])
                y.orig <- y
                if (missing(bwf)) {
                        N <- length(!is.na(y))
                        bw <- 1.06 * min(sd(y), IQR(y)/1.34) * 
                                N^-0.2
                }
                else {
                        bw <- bwf(y)
                        if (bw == 0) 
                                bw <- 0.1
                }
                if (max(y) <= 1 && 1 - max(y) < 2 * bw) {
                        if (min(y) >= 0 && min(y) < 2 * bw) {
                                scale <- "proportion"
                                y <- c(y, -y, 2 - y)
                        }
                }
                else if (min(y) >= 0 && min(y) < 2 * bw) {
                        scale <- "positive"
                        y <- c(y, -y)
                }
                else scale <- "open"
                dens <- density(y, bw = bw)
                if (scale == "proportion") {
                        dens$y[dens$x >= 0 & dens$x <= 1] <- 3 * 
                                dens$y[dens$x >= 0 & dens$x <= 
                                 1]
                        dens$y[dens$x < 0] <- 0
                        dens$y[dens$x > 1] <- 0
                        dens$x[dens$x < 0] <- 0
                        dens$x[dens$x > 1] <- 1
                }
                else if (scale == "positive") {
                        dens$y[dens$x >= 0] <- 2 * dens$y[dens$x >= 
                                0]
                        dens$y[dens$x < 0] <- 0
                        dens$x[dens$x < 0] <- 0
                }
                plot(dens, xlab = varnames(x)[j], ylab = "")
                points(y.orig, rep(0, length(y.orig)))
        }
        if (combine.chains) 
                for (j in 1:nvar(x)) {
                        do.densplot(j)
                        mtext(paste("Density of", varnames(x)[j]), 
                                side = 3, line = 2.5)
                        mtext(paste("(", format(niter(x)), " values per trace)", 
                                sep = ""), side = 3, line = 0.65)
                }
        else for (i in 1:nchain(x)) {
                for (j in 1:nvar(x)) {
                        do.densplot(j, i)
                        mtext(paste("Density from", chanames(x)[i]), 
                                side = 3, line = 2.5)
                        mtext(paste("(", format(niter(x)), " values per trace)", 
                                sep = ""), side = 3, line = 0.65)
                }
        }
        return()
}
"end.mcmc" <-
function (x) 
{
        attr(as.mcmc(x), "tspar")[2]
}
"frequency.mcmc" <-
function (x) 
{
        1/attr(as.mcmc(x), "tspar")[3]
}
"gelman.diag" <-
function (x, confidence = 0.95, transform = F) 
{
        #
        # Gelman and Rubin's code
        #
        # Adapted to work on mcmc objects. Now you can analyse
        # several variables at once.
        #
        x <- as.mcmc(x)
        #
        # We compute the following statistics:
        #
        #  xdot:  vector of sequence means
        #  s2:  vector of sequence sample variances (dividing by n-1)
        #  W = mean(s2):  within MS
        #  B = n*var(xdot):  between MS.
        #  muhat = mean(xdot):  grand mean; unbiased under strong stationarity
        #  varW = var(s2)/m:  estimated sampling var of W
        #  varB = B^2 * 2/(m-1):  estimated sampling var of B
        #  covWB = (n/m)*(cov(s2,xdot^2) - 2*muhat*cov(s^2,xdot)):
        #          estimated sampling cov(W,B)
        #  sig2hat = ((n-1)/n))*W + (1/n)*B:  estimate of sig2; unbiased under
        #            strong stationarity
        #  quantiles:  emipirical quantiles from last half of simulated
        #              sequences
        #
        if (nchain(x) == 1) 
                stop("You need at least two chains")
        if (start(x) < end(x)/2) 
                x <- window(x, start = end(x)/2 + 1)
        Niter <- niter(x)
        Nchain <- nchain(x)
        confshrink <- matrix(nrow = nvar(x), ncol = 2, dimnames = list(varnames(x), 
                c("Point est.", paste(50 * (1 + confidence), 
                        "% quantile", sep = ""))))
        for (i in 1:nvar(x)) {
                z <- x[, i, , drop = T]
                if (transform) 
                        if (min(z) > 0) 
                                z <- if (max(z) < 1) 
                                 log(z/(1 - z))
                                else log(z)
                s2 <- apply(z, 2, cov)
                W <- mean(s2)
                zbar <- apply(z, 2, mean)
                B <- Niter * var(zbar)
                sig2hat <- ((Niter - 1) * W + B)/Niter
                muhat <- mean(zbar)
                varW <- var(s2)/Nchain
                varB <- 2 * B^2/(Nchain - 1)
                covWB <- (Niter/Nchain) * (cov(s2, zbar^2) - 
                        2 * muhat * cov(s2, zbar))
                #
                # Posterior interval post.range combines all uncertainties
                # in a t interval with center muhat, scale sqrt(postvar),
                # and postvar.df degrees of freedom.
                #
                # postvar = sig2hat + B/(mn):  variance for the posterior
                #           interval. The B/(mn) term is there because of the
                #           sampling variance of muhat.
                # varpostvar:  estimated sampling variance of postvar
                #
                # 
                # Posterior interval post.range combines all uncertainties
                # in a t interval with center muhat, scale sqrt(postvar), 
                # and postvar.df degrees of freedom.
                #
                # postvar = sig2hat + B/(mn):  variance for the posterior
                #           interval. The B/(mn) term is there because of the
                #           sampling variance of muhat.
                # varpostvar:  estimated sampling variance of postvar
                #
                postvar <- sig2hat + B/(Niter * Nchain)
                varpostvar <- ((Niter - 1)^2 * varW + (1 + 1/Nchain)^2 * 
                        varB + 2 * (Niter - 1) * (1 + 1/Nchain) * 
                        covWB)/Niter^2
                post.df <- 2 * postvar^2/varpostvar
                df.adj <- (post.df + 3)/(post.df + 1)
                #
                # Estimated potential scale reduction (that would be achieved
                # by continuing simulations forever) has two components: an
                # estimate and an approx upper bound.
                #
                # confshrink = sqrt(postvar/W), 
                #     multiplied by sqrt((df+3)/(df+1)) as an adjustment for the
                #     width of the t-interval with df degrees of freedom.
                #
                # postvar/W = (n-1)/n + (1+1/m)(1/n)(B/W); we approximate
                # the sampling dist.  of (B/W) by an F distribution, with
                # degrees of freedom estimated from the approximate
                # chi-squared sampling dists for B and W.  (The F
                # approximation assumes that the sampling dists of B and W
                # are independent; if they are positively correlated, the
                # approximation is conservative.)
                #
                varlo.df <- 2 * W^2/varW
                R2.fixed <- (Niter - 1)/Niter
                R2.random <- (1 + 1/Nchain) * (1/Niter) * (B/W)
                R2.estimate <- R2.fixed + R2.random
                R2.upper <- R2.fixed + qf((1 + confidence)/2, 
                        Nchain - 1, varlo.df) * R2.random
                confshrink[i, 1] <- sqrt(df.adj * R2.estimate)
                confshrink[i, 2] <- sqrt(df.adj * R2.upper)
        }
        out <- list(confshrink = confshrink)
        class(out) <- "gelman.diag"
        out
}
"gelman.shrink.plot" <-
function (x, confidence = 0.95, transform = F, bin.width, max.bins, 
        auto.layout = T) 
{
        oldpar <- NULL
        if (auto.layout) 
                oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 
                        Nparms = nvar(x)))
        oldpar <- c(oldpar, par(set.scale()))
        on.exit(oldpar)
        if (!missing(bin.width)) {
                if (!missing(max.bins)) 
                        nbin <- min(floor((niter(x) - 50)/bin.width), 
                                max.bins)
                else nbin <- floor((niter(x) - 50)/bin.width)
        }
        else {
                if (!missing(max.bins)) 
                        nbin <- max.bins
                else nbin <- min(floor((niter(x) - 50)/10), 50)
        }
        bin <- floor((niter(x) - 50)/nbin)
        shrink <- array(dim = c(nrow = nbin + 1, ncol(x), 2))
        for (i in 1:nbin) {
                shrink[i, , ] <- gelman.diag(x[1:(50 + (i - 1) * 
                        bin), , ])$confshrink
        }
        shrink[nbin + 1, , ] <- gelman.diag(x)$confshrink
        title.scale <- par("cex") * switch(coda.options("ps.plot"), 
                1, 0.85, 0.7)
        scale <- par("cex") * switch(coda.options("ps.plot"), 
                1, 0.8, 0.6)
        all.na <- apply(is.na(shrink[, , 1, drop = F]), 2, all)
        if (any(all.na)) {
                cat("\n******* Error: *******\nCannot compute Gelman & Rubin's diagnostic for any chain \nsegments for variables", 
                        colnames(x)[all.na], "\nThis indicates convergence failure ==> Run chains for more iterations\n")
        }
        else for (j in 1:ncol(x)) {
                ymin <- min(c(1, shrink[, j, ]), na.rm = T)
                ymax <- max(c(1, shrink[, j, ]), na.rm = T)
                ylim <- range(pretty(c(ymin, ymax)))
                time.x <- unclass(time(x))[c(50 + 0:(nbin - 1) * 
                        bin, niter(x))]
                xlim <- range(pretty(time.x))
                plot(time.x, shrink[, j, 1], type = "l", xlab = "Last iteration in segment", 
                        ylab = "Shrink factor", cex = scale, 
                        xlim = xlim, ylim = ylim, axes = F)
                axis(1, at = pretty(c(50, 50 + (1:nbin) * bin) + 
                        i, 3), cex = scale, labels = format(pretty(c(50, 
                        50 + (1:nbin) * bin) + i, 3)))
                par(crt = 90, srt = 90)
                axis(2, at = pretty(c(ymin, ymax), 2), cex = scale, 
                        labels = format(pretty(c(ymin, ymax), 
                                2)))
                par(crt = 0, srt = 0)
                box()
                lines(time.x, shrink[, j, 2], lty = 2, col = "blue")
                abline(h = 1, lty = 4)
                i <- start(x) - 2
                legend(xlim[2], ylim[2], legend = c("median", 
                        paste(50 * (confidence + 1), "%", sep = "")), 
                        lty = c(1, 2), cex = scale * 0.7, bty = "n", 
                        col = c("black", "blue"), xjust = 1, 
                        yjust = 1)
                title(main = colnames(x)[j], cex = title.scale)
                if ((nna <- sum(is.na(shrink[, j, 1])))) {
                        cat("\n******* Warning: *******\nCould not compute Gelman & Rubin's diagnostic for", 
                                nna, "chain segments\n")
                }
                if (j != ncol(x) && mpause2()) 
                        break
        }
}
"gelman.trace.plot" <-
function (x, auto.layout = T) 
{
        # 
        # gandr.trace.plot -- function for plotting multiple chains on
        # same graph and printing Gelman & Rubin's convergence diagnostic
        # for each variable.
        # 
        # Author: Nicky Best 
        # 
        oldpar <- NULL
        if (auto.layout) 
                oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 
                        Nparms = nvar(x)))
        oldpar <- c(oldpar, par(set.scale()))
        on.exit(par(oldpar))
        title.scale <- par("cex") * switch(coda.options("ps.plot"), 
                0.8, 0.7, 0.65)
        scale <- par("cex") * switch(coda.options("ps.plot"), 
                1, 0.8, 0.6)
        for (i in 1:nvar(x)) {
                r <- x[, i, , drop = F]
                lim <- range(pretty(c(min(r), max(r))))
                plot.default(time(x), r[, , 1], type = "l", ylim = lim, 
                        xlab = "Iteration", ylab = colnames(r), 
                        cex = scale, axes = F)
                axis(1, at = pretty(time(x), 3), cex = scale, 
                        labels = format(pretty(time(x), 3)))
                par(crt = 90, srt = 90)
                axis(2, at = pretty(c(min(r), max(r)), 2), cex = scale, 
                        labels = format(pretty(c(min(r), max(r)), 
                                2)))
                par(crt = 0, srt = 0)
                box()
                for (j in 2:(dim(x)[3])) {
                        lines(r[, , j], lty = j, col = j)
                }
                shrink <- gelman.diag(r)$confshrink
                title(main = paste("Median =", format(round(shrink[1], 
                        2)), ",", "97.5% =", format(round(shrink[2], 
                        2))), cex = title.scale)
                if (i != nvar(x) && mpause2()) 
                        break
        }
}
"geweke.diag" <-
function (x, frac1 = 0.1, frac2 = 0.5) 
{
        xstart <- c(start(x), end(x) - frac2 * (end(x) - start(x)))
        xend <- c(start(x) + frac1 * (end(x) - start(x)), end(x))
        y.vom <- y.mean <- array(dim = c(2, nvar(x), nchain(x)), 
                dimnames = list(NULL, varnames(x), chanames(x)))
        for (i in 1:2) {
                y <- window(x, start = xstart[i], end = xend[i])
                y.mean[i, , ] <- apply(y, c(2, 3), mean)
                spans <- min(sqrt(niter(y))/0.3 + 1, niter(y) - 
                        1)
                for (j in 1:nvar(x)) for (k in 1:nchain(x)) {
                        spec0 <- spec.pgram(y[, j, k], spans = spans, 
                                demean = T, detrend = F, plot = F)$spec[1]
                        y.vom[i, j, k] <- 10^(spec0/10)/niter(y)
                }
        }
        z <- (y.mean[1, , ] - y.mean[2, , ])/sqrt(y.vom[1, , 
                ] + y.vom[2, , ])
        out <- list(z = z, frac = c(frac1, frac2))
        class(out) <- "geweke.diag"
        return(out)
}
"geweke.nse" <-
function (x) 
{
        # 
        # geweke.nse 
        # 
        # Author: Kate Cowles 
        # 
        # spans parm for smoothing periodogram 
        a <- length(x)
        #  
        nspans <- sqrt(a)/0.3 + 1
        if (nspans >= a) {
                # spans is longer than time series 
                nspans <- a - 1
        }
        pgram1 <- spec.pgram(x, spans = nspans, demean = T, detrend = F, 
                plot = F)
        # changed 03/20/94 
        # power <- 2 * pi * (10^(pgram1$spec[1]/10))     
        power <- (10^(pgram1$spec[1]/10))
        # spectral density converted from decibels       
        nse <- sqrt(power/a)
        nse
}
"geweke.plot" <-
function (x, frac1 = 0.1, frac2 = 0.5, max.bins, bin.width, auto.layout = T) 
{
        if (!missing(bin.width)) {
                if (!missing(max.bins)) 
                        nbin <- min(floor((niter(x) - 50)/bin.width), 
                                max.bins)
                else nbin <- floor((niter(x) - 50)/bin.width)
        }
        else {
                if (!missing(max.bins)) 
                        nbin <- max.bins
                else nbin <- min(floor((niter(x) - 50)/10), 50)
        }
        bin <- floor((niter(x) - 50)/nbin)
        oldpar <- NULL
        if (auto.layout) 
                oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 
                        Nparms = nvar(x)))
        oldpar <- c(oldpar, set.scale())
        on.exit(par(oldpar))
        title.scale <- par("cex") * c(1, 0.85, 0.7)[coda.options("ps.plot")]
        scale <- par("cex") * c(1, 0.8, 0.6)[coda.options("ps.plot")]
        ystart <- seq(from = start(x), to = end(x) - 49 * thin(x), 
                by = bin)
        gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 
                dimnames = c(ystart, varnames(x), chanames(x)))
        for (n in 1:length(ystart)) gcd[n, , ] <- geweke.diag(window(x, 
                start = ystart[n]), frac1 = frac1, frac2 = frac2)$z
        for (k in 1:nchain(x)) for (j in 1:nvar(x)) {
                ymax <- max(max(gcd[, j, k]), max(-gcd[, j, k]), 
                        2) * 1.4
                ymin <- -ymax
                plot(ystart, gcd[, j, k], type = "p", xlab = "First iteration in segment", 
                        ylab = "Z-score", cex = scale, ylim = c(ymin, 
                                ymax), pch = 4, axes = F)
                axis(1, at = pretty(ystart, 3), cex = scale, 
                        labels = format(pretty(ystart, 3)))
                par(crt = 90, srt = 90)
                axis(2, at = pretty(c(ymin, ymax), 2), cex = scale, 
                        labels = format(pretty(c(ymin, ymax), 
                                2)))
                par(crt = 0, srt = 0)
                box()
                abline(h = 1.96, lty = 2)
                abline(h = -1.96, lty = 2)
                if (dim(x)[3] > 1) {
                        title(main = paste(varnames(x)[j], " (", 
                                chanames(x)[k], ")", sep = ""), 
                                cex = title.scale)
                }
                else {
                        title(main = paste(varnames(x)[j], sep = ""), 
                                cex = title.scale)
                }
                if (!(k == nchain(x) && j == nvar(x)) && mpause2()) 
                        break
        }
        invisible(list(start.iter = ystart, z = gcd))
}
"geweke.power" <-
function (x) 
{
        # 
        # geweke.power 
        # 
        # spans parm for smoothing periodogram 
        a <- length(x)
        nspans <- sqrt(a)/0.3 + 1
        if (nspans >= a) {
                # spans is longer than time series 
                nspans <- a - 1
        }
        pgram1 <- spec.pgram(x, spans = nspans, demean = T, detrend = F, 
                plot = F)
        #changed 03/20/94        
        # power <- 2 * pi * (10^(pgram1$spec[1]/10))     
        power <- (10^(pgram1$spec[1]/10))
        # spectral density converted from decibels       
        nse <- sqrt(power/a)
        power
}
"heidel.diag" <-
function (x, eps = 0.1) 
{
        # 
        # Heidelberger and Welch diagnostic 
        # 
        # Author: Kate Cowles 
        # Modified by: Nicky Best 
        # 
        HW.list <- vector("list", nchain(x))
        names(HW.list) <- chanames(x)
        HW.mat0 <- matrix(0, ncol = 7, nrow = nvar(x))
        dimnames(HW.mat0) <- list(varnames(x), c("stest", "keep", 
                "discard", "C-vonM", "htest", "mean", "halfwidth"))
        for (i in 1:nchain(x)) {
                HW.mat <- HW.mat0
                for (j in 1:nvar(x)) {
                        Y <- x[, j, i, drop = T]
                        n1 <- length(Y)
                        n <- length(Y)
                        S0 <- geweke.power(Y[(n/2 + 1):n])
                        passed <- F
                        while (n >= n1/2 && !passed) {
                                T1 <- cumsum(Y)
                                ybar <- mean(Y)
                                B <- T1 - ybar * (1:n)
                                Bsq <- (B * B)/(n * S0)
                                I <- (2 * sum(Bsq[seq(2, n - 
                                 2, by = 2)]) + 4 * sum(Bsq[seq(1, 
                                 n - 1, by = 2)]) + Bsq[n])/(3 * 
                                 n)
                                passed <- !is.na(I) & I < 0.46
                                if (!passed) {
                                 Y <- Y[(n1/10 + 1):n]
                                 n <- length(Y)
                                }
                        }
                        S0ci <- geweke.power(Y)
                        halfwidth <- 1.96 * sqrt(S0ci/n)
                        passed2 <- (!is.na(halfwidth) & abs(halfwidth/ybar) <= 
                                eps)
                        if (is.na(I) | is.na(halfwidth) | !passed) {
                                n <- NA
                                nd <- NA
                                passed2 <- NA
                                ybar <- NA
                                halfwidth <- NA
                        }
                        else {
                                nd <- length(x[, j, i, drop = T]) - 
                                 n
                        }
                        HW.mat[j, ] <- c(passed, n * thin(x), 
                                nd * thin(x), I, passed2, ybar, 
                                halfwidth)
                }
                HW.list[[i]] <- HW.mat
        }
        class(HW.list) <- "heidel.diag"
        return(HW.list)
}
"[.mcmc" <-
function (x, i, j, k, drop = F) 
{
        #Drop = F by default, which is the opposite of the normal default
        #
        y <- NextMethod("[", drop = drop)
        if (drop) {
                attr(y, "tspar") <- NULL
                attr(y, "transform") <- NULL
                return(y)
        }
        else {
		if (!missing(j)) attr(y, "transform") <- attr(x, "transform")[j]
		else attr(y,"transform") <- attr(x,"transform")
		attr(y,"title") <- attr(x,"title")
	}
        if (missing(i)) 
                mcmc(y, start = start(x), thin = thin(x))
        else {
                n <- dim(x)[1]
                ind <- (1:n)[i]
                delta <- unique(ind[-1] - ind[-length(ind)])
                if (length(delta) != 1 || delta <= 0) {
                        warning("Not returning an MCMC object")
                }
                else {
                        xtspar <- tspar(x)
                        xtimes <- seq(from = xtspar[1], to = xtspar[2], 
                                by = xtspar[3])
                        ytspar <- numeric(3)
                        ytspar[1] <- xtimes[min(ind)]
                        ytspar[2] <- xtimes[max(ind)]
                        ytspar[3] <- (ytspar[2] - ytspar[1])/(length(ind) - 
                                1)
                        tspar(y) <- ytspar
                        class(y) <- "mcmc"
                }
                y
        }
}
"is.mcmc" <-
function (x) 
inherits(x, "mcmc")
"join.mcmc" <-
function (...) 
{
        x <- list(...)
        if (length(x) == 1 && is.list(x[[1]])) 
                x <- x[[1]]
        if (!all(unlist(lapply(x, is.mcmc)))) 
                stop("Arguments must be mcmc objects")
        nargs <- length(x)
        if (nargs == 1) 
                return(x[[1]])
        xtspar <- lapply(x, tspar)
        xtspar1 <- xtspar[[1]]
        if (!all(unlist(lapply(xtspar, "==", xtspar1)))) 
                stop("Arguments have different tspar attributes")
        xdim <- matrix(unlist(lapply(x, dim)), nrow = 3)
        if (any(xdim[1, ] != xdim[1, 1])) 
                stop("Chains have different numbers of iterations")
        if (any(xdim[2, ] != xdim[2, 1])) 
                stop("Chains have different numbers of variables")
        end.chn <- cumsum(xdim[3, ])
        start.chn <- c(1, end.chn[-nargs] + 1)
        n.chn <- sum(xdim[3, ])
        y <- array(NA, dim = c(xdim[1:2, 1], n.chn))
        ydn <- list(NULL, NULL, character(n.chn))
        for (i in 1:nargs) {
                y[, , start.chn[i]:end.chn[i]] <- x[[i]]
                xdn <- dimnames(x[[i]])
                if (!is.null(xdn)) {
                        if (!is.null(xdn[[1]])) 
                                if (is.null(ydn[[1]])) 
                                 ydn[[1]] <- xdn[[1]]
                                else if (any(ydn[[1]] != xdn[[1]])) 
                                 stop("Unequal variable names")
                        if (!is.null(xdn[[2]])) 
                                if (is.null(ydn[[2]])) 
                                 ydn[[2]] <- xdn[[2]]
                                else if (any(ydn[[2]] != xdn[[2]])) 
                                 stop("Unequal variable names")
                        if (!is.null(xdn[[3]])) 
                                ydn[[3]][start.chn[i]:end.chn[i]] <- xdn[[3]]
                        else ydn[[3]][start.chn[i]:end.chn[i]] <- rep(NA, 
                                dim(x[[i]])[3])
                }
        }
        dimnames(y) <- ydn
        mcmc(y, start = xtspar1[1], end = xtspar1[2], thin = xtspar1[3])
}
"mcmc" <-
function (data = NA, start = 1, end = numeric(0), thin = 1, ts.eps = .Options$ts.eps) 
{
        if (is.vector(data)) {
                ndata <- length(data)
                nvar <- 1
                nchains <- 1
                data <- array(data, c(ndata, nvar, nchains))
        }
        else if (is.matrix(data)) {
                ndata <- nrow(data)
                nvar <- ncol(data)
                nchains <- 1
                data.dimnames <- dimnames(data)
                data <- array(data, c(ndata, nvar, nchains))
                dimnames(data) <- c(data.dimnames, NA)
        }
        else if (is.array(data) && length(dim(data)) == 3) {
                ndata <- dim(data)[1]
                nvar <- dim(data)[2]
                nchains <- dim(data)[3]
        }
        else stop("Can\'t coerce data to mcmc object")
        if (length(start) > 1) 
                stop("Invalid start")
        if (length(end) > 1) 
                stop("Invalid end")
        if (length(thin) != 1) 
                stop("Invalid thin")
        if (missing(end)) 
                end <- start + (ndata - 1) * thin
        else if (missing(start)) 
                start <- end - (ndata - 1) * thin
        if (thin > 1 && abs(thin - round(thin)) < ts.eps) 
                thin <- round(thin)
        nobs <- floor((end - start)/thin + 1.01)
        if (ndata < nobs) 
                stop("Start, end and thin incompatible with data")
        else {
                end <- start + thin * (nobs - 1)
                if (nobs < ndata) 
                        data <- data[1:nobs, , , drop = FALSE]
        }
        attr(data, "tspar") <- c(start, end, thin)
        attr(data, "class") <- "mcmc"
        data
}
#Extractor functions for mcmc objects
#These render code a bit more legible and provide a degree of
#data abstraction

"niter" <- function(x) {
	if(is.mcmc(x))
		nrow(x)
}

"nvar" <- function(x) {
	if(is.mcmc(x))
		ncol(x)
}

"nchain" <- function(x) {
	if(is.mcmc(x))
		dim(x)[3]
}

"varnames" <- function(x) {
	if(is.mcmc(x))
		dimnames(x)[[2]]
}

"chanames" <- function(x) {
	if(is.mcmc(x))
		dimnames(x)[[3]]
}
"mpause2" <-
function () 
{
        mfg <- par("mfg")
        if (mfg[1] != mfg[3] || mfg[2] != mfg[4]) 
                return(F)
        cat("Save plot (s), quit plot (q) or continue (any other key):\n")
        ans <- readline()
        if (pmatch(ans, c("s", "S"), nomatch = 0)) {
                cat("filename:\n")
                filename <- readline()
                save.plot(filename)
                cat("Quit (q) or continue (any other key):\n")
                ans <- readline()
        }
        if (pmatch(ans, c("q", "Q"), nomatch = 0)) 
                return(T)
        else return(F)
}
"mpause3" <-
function () 
{
        cat("Save plot (s), or continue (any other key):\n")
        ans <- readline()
        if (pmatch(ans, c("s", "S"), nomatch = 0)) {
                cat("filename:\n")
                filename <- readline()
                save.plot(filename)
        }
}
"multi.menu" <-
function (choices, title, header, allow.zero = T) 
{
        # Select more than one value from a menu 
        # 
        if (!missing(title)) 
                cat(title, "\n\n")
        mat <- matrix(c(1:length(choices), choices), ncol = 2)
        if (!missing(header)) {
                if (length(header) == 2) 
                        mat <- rbind(header, mat)
                else stop("header is wrong length")
        }
        cat(paste(format(mat[, 1]), format(mat[, 2])), sep = "\n")
        repeat {
                cat("\nEnter relevant number(s), separated by commas", 
                        "Ranges such as 3:7 may be specified)", 
                        sep = "\n")
                if (allow.zero) 
                        cat("(Enter 0 for none)\n")
                ans <- scan(what = character(), sep = ",", strip.white = T, 
                        nlines = 1, quiet = T)
                if (length(ans) > 0) {
                        out <- numeric(0)
                        for (i in 1:length(ans)) {
                                nc <- nchar(ans[i])
                                wrd <- substring(ans[i], 1:nc, 
                                 1:nc)
                                colons <- wrd == ":"
                                err <- any(is.na(as.numeric(wrd[!colons]))) | 
                                 sum(colons) > 1 | colons[1] | 
                                 colons[nc]
                                if (err) {
                                 cat("Error: you have specified a non-numeric value!\n")
                                 break
                                }
                                else {
                                 out <- c(out, eval(parse(text = ans[i])))
                                 if (min(out) < ifelse(allow.zero, 
                                  0, 1) | max(out) > length(choices) | 
                                  (any(out == 0) & length(out) > 
                                   1)) {
                                  err <- T
                                  cat("Error: you have specified variable number(s) out of range!\n")
                                  break
                                 }
                                }
                        }
                        if (!err) 
                                break
                }
        }
        return(out)
}
"plot.mcmc" <-
function (x, trace = T, density = T, combine.chains = F, smooth = T, 
        bwf, auto.layout = T, one.page = F, ...) 
{
        oldpar <- NULL
        if (auto.layout) {
                mfrow <- set.mfrow(Nchains = nchain(x), Nparms = nvar(x), 
                        nplots = trace + density, sepplot = !combine.chains, 
                        one.page = one.page)
                oldpar <- par(mfrow = mfrow)
        }
        oldpar <- c(oldpar, par(set.scale()))
        on.exit(par(oldpar))
        if (combine.chains) {
                for (i in 1:nvar(x)) {
                        if (trace) 
                                traceplot(x[, i, , drop = F], 
                                 combine.chains = T, smooth = smooth)
                        if (density) 
                                if (missing(bwf)) 
                                 densplot(x[, i, , drop = F], 
                                  combine.chains = T)
                                else densplot(x[, i, , drop = F], 
                                 combine.chains = T, bwf = bwf)
                        if (i != nvar(x) && mpause2()) 
                                break
                }
        }
        else for (j in 1:nchain(x)) {
                for (i in 1:nvar(x)) {
                        if (trace) 
                                traceplot(x[, i, j], combine.chains = F, 
                                 smooth = smooth)
                        if (density) 
                                if (missing(bwf)) 
                                 densplot(x[, i, j, drop = F], 
                                  combine.chains = F)
                                else densplot(x[, i, j, drop = F], 
                                 combine.chains = F, bwf = bwf)
                        if (!(i == nvar(x) && j == nchain(x)) && 
                                mpause2()) {
                                stop.plot <- TRUE
                                break
                        }
                }
                if (exists("stop.plot", inherits = F) && stop.plot) 
                        break
        }
}
"plot.to.file" <-
function () 
{
        while (T) {
                cat("\nDo you want to save current plots as a postscript file (y/n) ?\n")
                ans <- scan(what = character(), nmax = 1, strip.white = T)
                if (length(ans) > 0) 
                        break
        }
        if (ans == "Y" | ans == "y") {
                while (T) {
                        cat("\nEnter name you want to call this postscript file:\n")
                        ps.name <- scan(what = character(), nmax = 1, 
                                strip.white = T)
                        if (length(ps.name) > 0) 
                                break
                }
                if (!is.null(Version()$language) && Version()$language == 
                        "R") 
                        save.plot(ps.name)
                else {
                        tmp <- switch(coda.options("ps.plot"), 
                                printgraph(method = "postscript", 
                                 file = ps.name, horizontal = T, 
                                 print = F, paper = "a4", width = 9, 
                                 height = 6, maximize = F), printgraph(method = "postscript", 
                                 file = ps.name, horizontal = F, 
                                 print = F, paper = "a4", width = 6, 
                                 height = 9, maximize = F), printgraph(method = "postscript", 
                                 file = ps.name, horizontal = F, 
                                 print = F, paper = "a4", width = 6, 
                                 height = 4, maximize = F))
                        if (is.null(tmp)) {
                                # default option 
                                printgraph(method = "postscript", 
                                 file = ps.name, horizontal = F, 
                                 print = F, paper = "a4", width = 6, 
                                 height = 9, maximize = F)
                        }
                }
        }
}
"print.coda.options" <-
function (data = F, stats = F, plots = F, diags = F) 
{
        #Display working data and coda options in pretty format 
        ## 
        # First define some formatting functions 
        # 
        doline <- function(x, title, type) {
                if (is.null(Version()$language)) 
                        doline <- get("doline", frame = sys.parent())
                if (length(x) > 1) {
                        doline(x[1], title)
                        for (i in 2:length(x)) doline(x[i], "")
                        return()
                }
                if (is.logical(x)) 
                        x <- ifelse(x, "Yes", "No")
                if (!missing(type) && mode(x) != type) 
                        x <- "N/A"
                len.title <- nchar(title)
                endchar <- ifelse(x == "", " ", ":")
                cat("| ", title, paste(rep(" ", 20 - nchar(title)), 
                        collapse = ""), sep = "")
                widthleft <- options("width")$width - 25 - nchar(x)
                if (widthleft > 0) 
                        cat("| ", x, paste(rep(" ", widthleft), 
                                collapse = ""), "|\n", sep = "")
                else cat("| ", x, "\n", sep = "")
        }
        dotitle <- function(title) {
                if (is.null(Version()$language)) 
                        doline <- get("doline", frame = sys.parent())
                doline("", title)
                doline("", paste(rep("-", nchar(title)), collapse = ""))
        }
        domaintitle <- function(title) {
                cat("| ", title, paste(rep(" ", options()$width - 
                        nchar(title) - 3), collapse = ""), "|\n", 
                        sep = "")
        }
        doblank <- function() {
                if (is.null(Version()$language)) 
                        doline <- get("doline", frame = sys.parent())
                doline("", "")
        }
        dosepline <- function() {
                cat("+", paste(rep("-", options("width")$width - 
                        2), collapse = ""), "+\n", sep = "")
        }
        strbrk <- function(x) {
                WIDTH <- options("width")$width - 26
                N <- length(x)
                x[-N] <- paste(x[-N], ", ", sep = "")
                tmp <- cumsum(nchar(x))
                lr <- vector("list", 1)
                lr[[1]] <- (1:N)[tmp <= WIDTH]
                rr <- x[lr[[1]]]
                i <- 1
                while (rr[length(rr)] != x[N]) {
                        tmp <- tmp + (WIDTH - sum(nchar(rr)))
                        lr <- c(lr, c(1:N)[tmp > (WIDTH * i) & 
                                tmp <= (WIDTH * (i + 1))])
                        rr <- x[lr[[i + 1]]]
                        i <- i + 1
                }
                y <- character(length(lr))
                for (i in 1:length(y)) y[[i]] <- paste(x[lr[[i]]], 
                        collapse = "")
                return(y)
        }
        # 
        # Now we can get on with it 
        # 
        cat("\nCurrent option settings:")
        cat("\n=======================\n\n")
        dosepline()
        if (data) {
                domaintitle("WORKING DATA")
                dosepline()
                doblank()
                ans <- strbrk(colnames(work.dat))
                doline(ans, "Variables selected")
                ans <- strbrk(dimnames(work.dat)[[3]])
                doline(ans, "Chains selected")
                wtspar <- tspar(work.dat)
                doline(wtspar[1], "Iterations - start")
                doline(wtspar[2], "             end")
                doline(wtspar[3], "Thinning interval")
                doblank()
                dosepline()
        }
        if (stats) {
                domaintitle("SUMMARY STATISTICS OPTIONS")
                dosepline()
                doblank()
                doline(coda.options("digits"), "Significant digits")
                doline(coda.options("combine.stats"), "Combine chains")
                doline(coda.options("batch.size"), "Batch size")
                ans <- strbrk(paste(coda.options("quantiles") * 
                        100, "%", sep = ""))
                doline(ans, "Quantiles")
                doblank()
                dosepline()
        }
        if (plots) {
                domaintitle("PLOTTING OPTIONS")
                dosepline()
                doblank()
                doline(coda.options("trace"), "Trace")
                doline(coda.options("lowess"), "Smooth lines")
                doline(coda.options("densplot"), "Density")
                width.cut <- options("width")[[1]] - 26
                if (is.null(options()$language)) {
                        func <- deparse(coda.options("bandwidth"))
                        for (i in 1:length(func)) {
                                func.i <- substring(func[i], 
                                 1:nchar(func[i]), 1:nchar(func[i]))
                                if (any(func.i == "\t")) {
                                 func.i[func.i == "\t"] <- "     "
                                 func[i] <- paste(func.i, collapse = "")
                                }
                        }
                }
                else func <- deparse(coda.options("bandwidth"), 
                        width.cutoff = width.cut)
                func.print <- vector("list", length(func))
                for (i in 1:length(func)) {
                        if (nchar(func[i]) <= width.cut) 
                                func.print[[i]] <- func[i]
                        else {
                                first <- seq(from = 1, to = nchar(func[i]), 
                                 by = width.cut)
                                last <- seq(from = width.cut, 
                                 to = nchar(func[i]), by = width.cut)
                                if (max(last) < nchar(func[i])) 
                                 last <- c(last, nchar(func[i]))
                                func.print[[i]] <- substring(func[i], 
                                 first = first, last = last)
                        }
                }
                doline(unlist(func.print), "Bandwidth")
                doline(coda.options("sepplot"), "Separate plot/chain")
                doline(coda.options("onepage"), "Separate page/var.")
                psopt <- switch(coda.options("ps.plot"), "Full page; Portrait", 
                        "Full page; Landscape", "Half page; Portrait")
                doline(psopt, "PostScript options")
                doblank()
                dosepline()
        }
        if (diags) {
                domaintitle("DIAGNOSTICS OPTIONS")
                dosepline()
                doblank()
                dotitle("Geweke")
                doline(coda.options("frac1"), "Window 1 fraction")
                doline(coda.options("frac2"), "Window 2 fraction")
                doline(coda.options("geweke.bin"), "Bin width", 
                        type = "numeric")
                doline(coda.options("geweke.max"), "Max number of bins", 
                        type = "numeric")
                doblank()
                dotitle("Gelman & Rubin")
                ans <- coda.options("gr.bin")
                doline(ans, "Bin width", type = "numeric")
                doline(coda.options("gr.max"), "Max number of bins", 
                        type = "numeric")
                doblank()
                dotitle("Raftery & Lewis")
                doline(coda.options("q"), "Quantile (q)")
                doline(coda.options("r"), "Precision (+/- r)")
                doline(coda.options("s"), "Probability (s)")
                doblank()
                dotitle("Cross-correlations")
                doline(coda.options("combine.corr"), "Combine chains")
                doblank()
                dosepline()
        }
        invisible()
}
"print.gelman.diag" <-
function (x, digits = 3, ...) 
{
        cat("Shrink factors:\n\n")
        print.default(x$confshrink, digits = digits, ...)
}
"print.geweke.diag" <-
function (x, digits = min(4, .Options$digits), ...) 
{
        cat("\nFraction in 1st window =", x$frac[1])
        cat("\nFraction in 2nd window =", x$frac[2], "\n\n")
        print.default(x$z, digits = digits, ...)
        invisible(x)
}
"print.heidel.diag" <-
function (x, digits = 3, ...) 
{
        HW.title <- matrix(c("Stationarity", "test", "# of iters.", 
                "to keep", "# of iters.", "to discard", "C-vonM", 
                "stat.", "Halfwidth", "test", "Mean", "", "Halfwidth", 
                ""), nrow = 2)
        for (i in 1:length(x)) {
                y <- matrix("", nrow = nrow(x[[i]]), ncol = 7)
                for (j in 1:ncol(y)) {
                        y[, j] <- format(x[[i]][, j], digits = digits, 
                                ...)
                }
                y[, c(1, 5)] <- ifelse(x[[i]][, c(1, 5)], "passed", 
                        "failed")
                cat("\n", names(x)[i], "\n\n", sep = "")
                y <- rbind(HW.title, y)
                dimnames(y) <- list(c("", "", dimnames(x[[i]])[[1]]), 
                        rep("", 7))
                print.default(y[, 1:4], quote = F, ...)
                print.default(y[, 5:7], quote = F, ...)
        }
        invisible(x)
}
"print.mcmc" <-
function (x, ...) 
{
        x <- as.mcmc(x)
        cat("Markov Chain Monte Carlo (MCMC) output:\nStart =", 
                start(x), "\nEnd =", end(x), "\nThinning interval =", 
                thin(x), "\nNumber of chains =", dim(x)[3], "\n\n")
        if (is.null(dimnames(x))) 
                dimnames(x) <- list(NULL, NULL, NULL)
        if (is.null(dimnames(x)[[1]])) 
                dimnames(x)[[1]] <- as.vector(time(x))
        if (is.null(dimnames(x)[[3]])) 
                dimnames(x)[[3]] <- paste("chain", 1:dim(x)[3])
        attr(x, "tspar") <- NULL
        attr(x, "class") <- NULL
        for (i in 1:dim(x)[3]) {
                cat(dimnames(x)[[3]][i], "\n\n", sep = "")
                print.default(array(x[, , i, drop = T], dim = dim(x)[1:2], 
                        dimnames = dimnames(x)[1:2]))
        }
        invisible(x)
}
"print.raftery.diag" <-
function (x, ...) 
{
        cat("\nRAFTERY AND LEWIS CONVERGENCE DIAGNOSTIC:")
        cat("\n=========================================\n\n")
        cat("Iterations used = ", x$tspar[1], ":", x$tspar[2], 
                "\n", sep = "")
        cat("Thinning interval =", x$tspar[3], "\n")
        cat("Sample size per chain =", x$Niters, "\n")
        cat("\nQuantile (q) =", x$params["q"])
        cat("\nAccuracy (r) = +/-", x$params["r"])
        cat("\nProbability (s) =", x$params["s"], "\n")
        if (x$resmatrix[1] == "Error") 
                cat("\nYou need a sample size of at least", x$resmatrix[2], 
                        "with these values of q, r and s\n")
        else for (i in 1:(dim(x$resmatrix)[3])) {
                cat("\n", dimnames(x$resmatrix)[[3]][i], "\n", 
                        sep = "")
                out <- x$resmatrix[, , i, drop = T]
                out <- rbind(matrix(c("Burn-in ", "Total", 
                        "Lower bound ", "Dependence", 
                        "(M)", "(N)", "(Nmin)", "factor (I)"), 
                        byrow = T, ncol = 4), out)
                dimnames(out) <- list(dimnames(out)[[1]], rep("", 
                        4))
                for (i in ncol(out)) out[, i] <- format(out[, 
                        i])
                print.default(out, quote = F)
        }
	invisible(x)
}
"print.summary.mcmc" <-
function (x, digits = max(3, .Options$digits - 3), ...) 
{
        cat("Iterations used = ", x$mcpar["start"], ":", x$mcpar["end"], 
                "\n", sep = "")
        cat("Thinning interval =", x$mcpar["thin"], "\n")
        cat("Sample size per chain =", (x$mcpar["end"] - x$mcpar["start"])/x$mcpar["thin"] + 
                1, "\n")
        if (x$combined) {
                cat("\nPooling over chains:\n", paste(x$chain.names, 
                        collapse = "\n"), "\n", sep = "")
        }
        cat("\n1. Empirical mean and standard deviation for each variable,")
        cat("\n   plus standard error of the mean:\n\n")
        if (length(dim(x$statistics)) == 3) 
                for (i in 1:(dim(x$statistics)[3])) {
                        cat("\n", x$chain.names[i], "\n", sep = "")
                        print(x$statistics[, , i], digits = digits, 
                                ...)
                }
        else print(x$statistics, digits = digits, ...)
        cat("\n2. Quantiles for each variable:\n\n")
        if (length(dim(x$quantiles)) == 3) 
                for (i in 1:(dim(x$quantiles)[3])) {
                        cat("\n", x$chain.names[i], "\n", sep = "")
                        print(x$quantiles[, , i], digits = digits, 
                                ...)
                }
        else print(x$quantiles, digits = digits, ...)
        invisible(x)
}
"raftery.diag" <-
function (data, q = 0.025, r = 0.005, s = 0.95, converge.eps = 0.001) 
{
        Nparms <- dim(data)[2]
        Niters <- dim(data)[1]
        Nchains <- dim(data)[3]
        resmatrix <- array(dim = c(Nparms, 4, Nchains))
        # 
        dimnames(resmatrix) <- list(dimnames(data)[[2]], c("M", 
                "N", "Nmin", "I"), dimnames(data)[[3]])
        phi <- qnorm(0.5 * (1 + s))
        nmin <- as.integer(ceiling((q * (1 - q) * phi^2)/r^2))
        if (nmin > Niters) 
                resmatrix <- c("Error", nmin)
        else for (k in 1:Nchains) {
                for (i in 1:Nparms) {
                        #          First need to find the thinning parameter kthin 
                        # 
                        quant <- quantile(data[, i, k, drop = T], 
                                probs = q)
                        dichot <- data[, i, k, drop = T] <= quant
                        kwork <- 0
                        bic <- 1
                        while (bic >= 0) {
                                kwork <- kwork + 1
                                testres <- dichot[seq(1, nrow(data), 
                                 by = kwork)]
                                newdim <- length(testres)
                                testtran <- table(testres[1:(newdim - 
                                 2)], testres[2:(newdim - 1)], 
                                 testres[3:newdim])
                                g2 <- 0
                                for (i1 in 1:2) {
                                 for (i2 in 1:2) {
                                  for (i3 in 1:2) {
                                   if (testtran[i1, i2, i3] != 
                                    0) {
                                    fitted <- (sum(testtran[i1, 
                                     i2, 1:2]) * sum(testtran[1:2, 
                                     i2, i3]))/(sum(testtran[1:2, 
                                     i2, 1:2]))
                                    g2 <- g2 + testtran[i1, i2, 
                                     i3] * log(testtran[i1, i2, 
                                     i3]/fitted) * 2
                                   }
                                  }
                                 }
                                }
                                bic <- g2 - log(newdim - 2) * 
                                 2
                        }
                        kthin <- as.integer(kwork * thin(data))
                        #
                        # then need to find length of burn-in and No of iterations for required precision 
                        # 
                        finaltran <- table(testres[1:(newdim - 
                                1)], testres[2:newdim])
                        alpha <- finaltran[1, 2]/(finaltran[1, 
                                1] + finaltran[1, 2])
                        beta <- finaltran[2, 1]/(finaltran[2, 
                                1] + finaltran[2, 2])
                        tempburn <- log((converge.eps* (alpha + beta))/max(alpha, 
                                beta))/(log(abs(1 - alpha - beta)))
                        nburn <- as.integer(ceiling(tempburn) * 
                                kthin)
                        tempprec <- ((2 - alpha - beta) * alpha * 
                                beta * phi^2)/(((alpha + beta)^3) * 
                                r^2)
                        nkeep <- as.integer(ceiling(tempprec) * 
                                kthin)
                        iratio <- (nburn + nkeep)/nmin
                        resmatrix[i, 1, k] <- nburn
                        resmatrix[i, 2, k] <- nkeep + nburn
                        resmatrix[i, 3, k] <- nmin
                        resmatrix[i, 4, k] <- signif(iratio, 
                                digits = 3)
                }
        }
        y <- list(tspar = tspar(data), params = c(r = r, s = s, 
                q = q), Niters = Niters, resmatrix = resmatrix)
        class(y) <- "raftery.diag"
        return(y)
}
"read.and.check" <-
function (message = "", what = numeric(), lower, upper, answer.in, 
        default) 
{
        #Read data from the command line and check that it satisfies 
        #certain conditions.  The function will loop until it gets 
        #and answer satisfying the conditions. This entails extensive 
        #checking of the conditions to  make sure they are consistent 
        #so we don't end up in an infinite loop. 
        have.lower <- !missing(lower)
        have.upper <- !missing(upper)
        have.ans.in <- !missing(answer.in)
        have.default <- !missing(default)
        if (have.lower | have.upper) {
                if (!is.numeric(what)) 
                        stop("Can't have upper or lower limits with non numeric input")
                if (have.lower && !is.numeric(lower)) 
                        stop("lower limit not numeric")
                if (have.upper && !is.numeric(upper)) 
                        stop("upper limit not numeric")
                if ((have.upper & have.lower) && upper < lower) 
                        stop("lower limit greater than upper limit")
        }
        if (have.ans.in) {
                if (mode(answer.in) != mode(what)) 
                        stop("inconsistent values of what and answer.in")
                if (have.lower) 
                        answer.in <- answer.in[answer.in >= lower]
                if (have.upper) 
                        answer.in <- answer.in[answer.in <= upper]
                if (length(answer.in) == 0) 
                        stop("No possible response matches conditions")
        }
        if (have.default) {
                if (mode(default) != mode(what)) 
                        stop("inconsistent values of what and default")
                if (have.lower && default < lower) 
                        stop("default value below lower limit")
                if (have.upper && default > upper) 
                        stop("default value above upper limit")
                if (have.ans.in && !any(answer.in == default)) 
                        stop("default value does not satisfy conditions")
        }
        err <- T
        while (err) {
                if (nchar(message) > 0) {
                        cat("\n", message, "\n", sep = "")
                        if (have.default) 
                                cat("(Default = ", default, ")\n", 
                                 sep = "")
                }
                else cat("\n")
                repeat {
                        cat("1:")
                        ans <- readline()
                        if (length(ans) == 1 && nchar(ans) > 
                                0) 
                                break
                        else if (have.default) {
                                ans <- default
                                break
                        }
                }
                if (is.numeric(what)) {
                        err1 <- T
                        ans <- as.numeric(ans)
                        message <- "You must enter a number"
                        if (is.na(ans)) 
                                NULL
                        else if ((have.lower & have.upper) && 
                                (ans < lower | ans > upper)) 
                                message <- paste(message, "between", 
                                 lower, "and", upper)
                        else if (have.lower && ans < lower) 
                                message <- paste(message, ">=", 
                                 lower)
                        else if (have.upper && ans > upper) 
                                message <- paste(message, "<=", 
                                 upper)
                        else err1 <- F
                }
                else err1 <- F
                if (have.ans.in) {
                        if (!is.na(ans) && !any(ans == answer.in)) {
                                message <- paste("You must enter one of the following:", 
                                 paste(answer.in, collapse = ","))
                                err2 <- T
                        }
                        else err2 <- F
                }
                else err2 <- F
                err <- err1 | err2
        }
        return(ans)
}
"read.bugs" <-
function (file = "bugs.out", start, end, thin) 
{
        #Based on readdat by Karen Vines 
        #Index object not required  
        #We don't remove missing values - each row of the output should 
        #contain at least one non-missing value. 
        #We return an mcmc object, if possible 
        # 
        nc <- nchar(file)
        if (nc > 3 && substring(file, nc - 3, nc) == ".out") 
                root <- substring(file, 1, nc - 4)
        else root <- file
        index <- read.table(file = paste(root, ".ind", sep = ""), 
                row.names = 1, col.names = c("", "begin", "end"))
        vnames <- row.names(index)
        # 
        temp <- scan(file = paste(root, ".out", sep = ""), what = list(iter = 0, 
                val = 0), quiet = T)
        # Do one pass through the data to see if we can construct 
        # a regular time series easily 
        # 
        start.vec <- end.vec <- thin.vec <- numeric(nrow(index))
        for (i in 1:length(vnames)) {
                iter.i <- temp$iter[index[i, "begin"]:index[i, 
                        "end"]]
                thin.i <- unique(diff(iter.i))
                thin.vec[i] <- if (length(thin.i) == 1) 
                        thin.i
                else NA
                start.vec[i] <- iter.i[1]
                end.vec[i] <- iter.i[length(iter.i)]
        }
        if (any(is.na(start.vec)) || any(thin.vec != thin.vec[1]) || 
                any((start.vec - start.vec[1])%%thin.vec[1] != 
                        0)) {
                # 
                # Do it the brute force way 
                # 
                iter <- sort(unique(temp$iter))
                old.thin <- unique(diff(iter))
                if (length(old.thin) == 1) 
                        is.regular <- TRUE
                else {
                        if (all(old.thin%%min(old.thin) == 0)) 
                                old.thin <- min(old.thin)
                        else old.thin <- 1
                        is.regular <- FALSE
                }
        }
        else {
                iter <- seq(from = min(start.vec), to = max(end.vec), 
                        by = thin.vec[1])
                old.thin <- thin.vec[1]
                is.regular <- TRUE
        }
        if (missing(start)) 
                start <- min(start.vec)
        else if (start < min(start.vec)) {
                warning("start not changed")
                start <- min(start.vec)
        }
        else if (start > max(end.vec)) 
                stop("Start after end of data")
        else iter <- iter[iter >= start]
        if (missing(end)) 
                end <- max(end.vec)
        else if (end > max(end.vec)) {
                warning("end not changed")
                end <- max(end.vec)
        }
        else if (end < min(start.vec)) 
                stop("End before start of data")
        else iter <- iter[iter <= end]
        if (missing(thin)) 
                thin <- old.thin
        else if (thin%%old.thin != 0) {
                thin <- old.thin
                warning("thin not changed")
        }
        else {
                new.iter <- iter[(iter - start)%%thin == 0]
                new.thin <- unique(diff(new.iter))
                if (length(new.thin) != 1 || new.thin != thin) 
                        warning("thin not changed")
                else {
                        iter <- new.iter
                        end <- max(iter)
                        is.regular <- TRUE
                }
        }
        out <- matrix(NA, nrow = length(iter), ncol = nrow(index))
        dimnames(out) <- list(iter, vnames)
        for (v in vnames) {
                cat("Abstracting", v, "... ")
                inset <- index[v, "begin"]:index[v, "end"]
                iter.v <- temp$iter[inset]
                if (!is.regular) {
                        use.v <- duplicated(c(iter, iter.v))[-(1:length(iter))]
                        use <- duplicated(c(iter.v, iter))[-(1:length(iter.v))]
                }
                else {
                        use.v <- (iter.v - start)%%thin == 0 & 
                                iter.v >= start & iter.v <= end
                        use <- (iter.v[use.v] - start)%/%thin + 
                                1
                }
                if (any(use) & any(use.v)) 
                        out[use, v] <- temp$val[inset[use.v]]
                cat(length(use), "valid values\n")
        }
        if (is.regular) {
                out <- mcmc(out, start = start, end = end, thin = thin)
                dimnames(out)[3] <- file
        }
        else warning("not returning an mcmc object")
        return(out)
}
"read.bugs.interactive" <-
function () 
{
        #Interactively read BUGS output 
        # 
        cat("\nEnter problem title, followed by return key:\n")
        mtitle <- scan(what = character(), sep = "\n", nmax = 1, 
                strip.white = T, quiet = T)
        cat("Enter BUGS output filenames, separated by return key\n")
        cat("(Leave blank and hit return to exit):\n")
        filenames <- scan(what = character(), sep = "\n", strip.white = T, 
                quiet = T)
        if (length(filenames) == 0) {
                warning("No file name given")
                cat("Using default name \"bugs\"\n")
                filenames <- "bugs"
        }
        nfiles <- length(filenames)
        chains <- vector("list", nfiles)
        names(chains) <- filenames
        tspar.mat <- matrix(NA, nrow = nfiles, ncol = 3)
        niters <- vector("numeric", length = nfiles)
        allnames <- character(0)
        for (i in 1:nfiles) {
                chains[[i]] <- read.bugs(file = filenames[i])
                tspar.mat[i, ] <- tspar(chains[[i]])
                niters[i] <- nrow(chains[[i]])
                allnames <- unique(c(allnames, colnames(chains[[i]])))
        }
        if (any(tspar.mat[, 2] != tspar.mat[1, 2])) 
                stop("Chains have different thinning intervals\nYou cannot combine them")
        if (nfiles > 1) {
                parm.OK <- rep(F, length(allnames))
                for (i in 1:length(allnames)) for (j in 1:nfiles) {
                        parm.OK[i] <- any(colnames(chains[[j]]) == 
                                allnames[i])
                        if (!parm.OK[i]) 
                                break
                }
                if (all(parm.OK == F)) {
                        cat("\n******* Warning: *******\nThere are no variables common to all chains\nOnly first \nchain (", 
                                filenames[1], ".out) will be used\n", 
                                sep = "")
                        chains <- chains[[1]]
                        filenames <- filenames[1]
                        nfiles <- 1
                }
                else if (any(parm.OK == F)) {
                        warning(paste("Only variable(s) ", paste(allnames[parm.OK], 
                                sep = ", "), " are common to all chains"))
                        cat("All others will be deleted\n")
                }
                for (i in 1:nfiles) {
                        chains[[i]] <- chains[[i]][, allnames[parm.OK], 
                                ]
                }
                if (!all(niters == niters[1])) {
                        warning("Chains are of different lengths\n")
                        cat("Truncating each to last", min(niters), 
                                "iterations\n")
                        for (i in 1:nfiles) {
                                chains[[i]] <- chains[[i]][(niters[i] - 
                                 min(niters) + 1):niters[i], 
                                 , ]
                        }
                }
                if (!all(tspar.mat[, 1] == tspar.mat[1, 1])) {
                        cat("\n******* Warning:  *******\nChains start at different iteration numbers\nAll plots will \nbe labelled by iteration numbers for first chain (", 
                                filenames[1], ")\n", sep = "")
                }
        }
        y <- join.mcmc(chains)
        attr(y, "title") <- mtitle
        return(y)
}
"set.mfrow" <-
function (Nchains = 1, Nparms = 1, nplots = 1, sepplot = F, one.page = F) 
{
        #  
        # Set up dimensions of graphics window: 
        # If only density plots OR trace plots are requested, dimensions are: 
        #	1 x 1	if Nparms = 1 
        #	1 X 2 	if Nparms = 2 
        #	2 X 2 	if Nparms = 3 or 4 
        #	3 X 2 	if Nparms = 5 or 6 or 10 - 12 
        #	3 X 3 	if Nparms = 7 - 9 or >= 13 
        # If both density plots AND trace plots are requested, dimensions are: 
        #	1 x 2	if Nparms = 1 
        #	2 X 2 	if Nparms = 2 
        #	3 X 2 	if Nparms = 3, 5, 6, 10, 11, or 12 
        #	4 x 2	if Nparms otherwise 
        # If separate plots are requested for each chain, dimensions are: 
        #	1 x 2	if Nparms = 1 & Nchains = 2 
        #	2 X 2 	if Nparms = 2 & Nchains = 2 OR Nparms = 1 & Nchains = 3 or 4 
        #	3 x 2	if Nparms = 3 or >= 5 & Nchains = 2  
        #		   OR Nchains = 5 or 6 or 10 - 12 (and any Nparms) 
        #	2 x 3	if Nparms = 2 or 4 & Nchains = 3 
        #	4 x 2   if Nparms = 4 & Nchains = 2  
        #		   OR Nchains = 4 & Nparms > 1 
        #	3 x 3	if Nparms = 3 or >= 5  & Nchains = 3  
        #		   OR Nchains = 7 - 9 or >= 13 (and any Nparms) 
        # If more plots are required than will fit on one page, the 
        # browser function is used to cycle round multiple pages of plots 
        #  
        if (one.page) 
                Nparms <- 1
        if (sepplot && Nchains > 1) {
                if (nplots == 1) {
                        mrows <- ifelse((Nparms == 1 && Nchains == 
                                2), 1, ifelse((Nparms == 2 && 
                                Nchains == 2) || (Nparms == 1 && 
                                any(Nchains == c(3, 4))) || (any(Nparms == 
                                c(2, 4)) && Nchains == 3), 2, 
                                ifelse((Nparms == 4 && Nchains == 
                                 2) || (Nchains == 4 && Nparms > 
                                 1), 4, 3)))
                        mcols <- ifelse((Nchains == 3 && any(Nparms == 
                                c(2, 4))) || ((Nparms == 3 || 
                                Nparms >= 5) && Nchains == 3) || 
                                any(Nchains == c(7:9)) || Nchains >= 
                                13, 3, 2)
                }
                else {
                        mrows <- ifelse(Nchains <= 4, Nchains, 
                                ifelse(any(Nchains == c(5, 6, 
                                 10:12)), 3, 4))
                        mcols <- 2
                }
        }
        else {
                if (nplots == 1) {
                        mrows <- ifelse(Nparms <= 2, 1, ifelse(Nparms <= 
                                4, 2, 3))
                        mcols <- ifelse(Nparms <= 9, Nparms%/%mrows + 
                                ifelse(Nparms%%mrows == 0, 0, 
                                 1), ifelse(any(Nparms == c(10:12)), 
                                2, 3))
                }
                else {
                        mrows <- ifelse(Nparms <= 4, Nparms, 
                                ifelse(any(Nparms == c(5, 6, 
                                 10:12)), 3, 4))
                        mcols <- 2
                }
        }
        return(mfrow = c(mrows, mcols))
}
"set.scale" <-
function () 
{
        set.scale.windows()
}
"set.scale.motif" <-
function (mrows, mcols) 
{
        plot.dim <- mrows + mcols - 1
        if (plot.dim > 7) 
                plot.dim <- 7
        if (mrows >= 6 & mcols < 6) 
                plot.dim <- 8
        if (mrows < 6 & mcols >= 6) 
                plot.dim <- 9
        if (mrows >= 6 & mcols >= 6) 
                plot.dim <- 10
        par(mfrow = c(mrows, mcols), oma = c(1.1, 0, 3.8, 0), 
                mgp = c(3.9, 0.9, 0))
        par(mar = c(6.1, 6.5, 5.1, 3.1))
        args.list$scale <- switch(plot.dim, 1.1, 1, 0.8, 0.75, 
                0.65, 0.6, 0.6, 0.35, 0.55, 0.35)
        args.list$title.scale <- switch(plot.dim, 0.825, 0.775, 
                0.7, 0.55, 0.5, 0.5, 0.5, 0.45, 0.3, 0.3)
        if (plot.dim < 3 & coda.options("ps.plot") < 3) {
                par(lab = c(4, 4, 5))
        }
        else {
                if (plot.dim < 6) {
                        par(lab = c(4, 2, 4))
                }
                else {
                        par(lab = c(3, 2, 3))
                }
        }
        return(args.list)
}
"set.scale.windows" <-
function () 
{
        mrows <- par("mfrow")[1]
        mcols <- par("mfrow")[2]
        plot.dim <- mrows + mcols - 1
        if (plot.dim > 7) 
                plot.dim <- 7
        if (mrows >= 6 & mcols < 6) 
                plot.dim <- 8
        if (mrows < 6 & mcols >= 6) 
                plot.dim <- 9
        if (mrows >= 6 & mcols >= 6) 
                plot.dim <- 10
        oma <- c(1, 0, 3.8, 0)
        mgp <- c(1.5, 0.4, 0)
        mar <- c(3.25, 3.5, 1.5, 2)
        cex <- switch(plot.dim, 0.9, 0.8, 0.7, 0.65, 0.6, 0.5, 
                0.425, 0.375, 0.2, 0.2)
        lab <- if (plot.dim < 3 & coda.options("ps.plot") < 3) {
                c(4, 4, 5)
        }
        else {
                if (plot.dim < 6) 
                        c(4, 2, 4)
                else c(3, 2, 3)
        }
        return(oma = oma, mgp = mgp, mar = mar, cex = cex, lab = lab)
}
"set.scale.x11" <-
function () 
{
        mrows <- par("mfrow")[1]
        mcols <- par("mfrow")[2]
        plot.dim <- mrows + mcols - 1
        if (plot.dim > 7) 
                plot.dim <- 7
        if (mrows >= 6 & mcols < 6) 
                plot.dim <- 8
        if (mrows < 6 & mcols >= 6) 
                plot.dim <- 9
        if (mrows >= 6 & mcols >= 6) 
                plot.dim <- 10
        oma <- c(1.3, 0, 4.2, 0)
        mgp <- c(4.2, 0.9, 0)
        mar <- c(6.7, 7, 5.7, 3.5)
        cex <- switch(plot.dim, 1.1, 1, 0.8, 0.75, 0.65, 0.6, 
                0.6, 0.35, 0.55, 0.35)
        lab <- if (plot.dim < 3 & coda.options("ps.plot") < 3) {
                c(4, 4, 5)
        }
        else {
                if (plot.dim < 6) 
                        c(4, 2, 4)
                else c(3, 2, 3)
        }
        return(oma = oma, mgp = mgp, mar = mar, cex = cex, lab = lab)
}
"spec.pgram" <-
function (x, spans = 1, taper = 0.1, demean = F, detrend = T, 
        pad = F, plot = F) 
{
        N <- length(x)
        if (detrend) {
                t <- 1:N
                x <- residuals(lm(x ~ t))
        }
        else if (demean) {
                x <- x - mean(x)
        }
        if (taper > 0.5 || taper < 0) 
                stop("taper must be between 0 and 0.5")
        else if (taper > 0) {
                w <- rep(1, N)
                n <- max(round(N * taper), 1)
                w[1:n] <- sin(((1:n - 0.5) * pi)/(2 * n))^2
                w[N:(N - n + 1)] <- w[1:n]
                x <- x * w
        }
        if (pad) 
                x <- c(x, rep(0, nextn(N) - N))
        Nspec <- ceiling(N/2) + 1
        spec <- (Mod(fft(x))^2)[1:Nspec]/N
        for (i in spans) {
                m <- floor(i/2)
                if (m > 0) {
                        filter <- c(0.5, rep(1, 2 * m - 1), 0.5)/(2 * 
                                m)
                        spec0 <- c(spec[(m + 1):2], spec, spec[Nspec - 
                                1:m])
                        for (j in 1:Nspec) spec[j] <- sum(spec0[j:(j + 
                                2 * m)] * filter)
                }
        }
        spec <- 10 * log10(spec)
        freq <- seq(from = 0, to = 0.5, length = Nspec)
        return(spec = spec, freq = freq)
}
"start.mcmc" <-
function (x) 
{
        attr(as.mcmc(x), "tspar")[1]
}
"summary.mcmc" <-
function (x, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), combine.chains = F, 
        ...) 
{
        chain.names <- chanames(x)
        mcpar <- c(start = start(x), end = end(x), thin = thin(x))
        combine.chains <- combine.chains && (nchain(x) > 1)
        if (combine.chains) 
                Nchains.stats <- 1
        else Nchains.stats <- nchain(x)
        statnames <- c("Mean", "SD", "Naive SE", "Time-series SE")
        chain.names <- if (combine.chains) 
                NULL
        else chanames(x)
        varstats <- array(dim = c(nvar(x), 4, Nchains.stats), 
                dimnames = list(varnames(x), statnames, chain.names))
        for (j in 1:Nchains.stats) {
                if (combine.chains) {
                        xmean <- apply(x, 2, mean, na.rm = T)
                        xvar <- apply(x, 2, var, na.rm = T)
                        xvar <- diag(var(matrix(aperm(x, c(2, 
                                1, 3)), ncol = nvar(x), byrow = T)))
                        xnse <- apply(x, 2, geweke.nse)
                        n.used <- niter(x) * nchain(x)
                }
                else {
                        xmean <- apply(x[, , j], 2, mean, na.rm = T)
                        xvar <- apply(x[, , j], 2, var, na.rm = T)
                        xnse <- apply(x[, , j], 2, geweke.nse)
                        n.used <- niter(x)
                }
                varstats[, 1, j] <- xmean
                varstats[, 2, j] <- sqrt(xvar)
                varstats[, 3, j] <- sqrt(xvar/n.used)
                varstats[, 4, j] <- xnse
        }
        varstats <- drop(varstats)
        varquant <- if (combine.chains) 
                t(apply(x, 2, quantile, quantiles))
        else aperm(apply(x, c(2, 3), quantile, quantiles), c(2, 
                1, 3))
        varquant <- drop(varquant)
        out <- list(statistics = varstats, quantiles = varquant, 
                mcpar = mcpar, combined = combine.chains, chain.names = chanames(x))
        class(out) <- "summary.mcmc"
        return(out)
}
"test.bandwidth" <-
function () 
{
        # 
        # test.bandwidth -- Attempts to evaluate the bandwidth function specified 
        #                   by the user. If a single numeric value results, then 
        #                   the function is passed, else an error is returned or 
        #                   the function is dumped  
        # 
        # Author: Nicky Best 
        # 
        out <- F
        y <- 1:10
        err <- coda.options("bandwidth")(y)
        if (is.numeric(err) & length(err) == 1) {
                out <- T
        }
        else {
                cat("Error: this function is not appropriate for calculating kernel bandwidths!\n\n")
                out <- F
        }
        out
}
"thin" <-
function (x, ...) 
UseMethod("thin")
"thin.mcmc" <-
function (x) 
{
        attr(as.mcmc(x), "tspar")[3]
}
"tidy.up" <-
function () 
{
        # 
        # tidy.up -- gives option to save input files in S format; then deletes all  
        #	     S-plus objects created during current CODA session, and 
        #	     closes all graphics windows 
        # 
        # Author: Nicky Best 
        # 
        cat("\nQuitting CODA....\n")
        if (exists("coda.dat", where = 1)) {
                while (T) {
                        cat("\nDo you want to save the BUGS output as an R object file(y/n) ?\n")
                        ans <- scan(what = character(), nmax = 1, 
                                strip.white = T)
                        if (length(ans) > 0) 
                                break
                }
                if (ans == "Y" | ans == "y") {
                        cat("Enter name you want to call this object file:\n")
                        fname <- scan(what = character(), nmax = 1, 
                                strip.white = T)
                        assign(fname, coda.dat, envir = sys.frame(0))
                }
        }
        coda.objects <- c("coda.dat", "work.dat")
        if (!is.null(Version()$language) && Version()$language == 
                "R") 
                for (i in coda.objects) {
                        if (exists(i)) 
                                rm(list = i, inherits = TRUE)
                }
        else if (is.null(Version()$language)) {
                for (i in coda.objects) {
                        if (exists(i)) 
                                remove(i, where = 1)
                }
                graphics.off()
        }
        else error("language not recognized")
        cat("Have a BUGS-free day!\n")
}
"time.mcmc" <-
function (x) 
{
        x <- as.mcmc(x)
        xtspar <- attr(x, "tspar")
        ts(seq(from = xtspar[1], to = xtspar[2], by = xtspar[3]), 
                start = xtspar[1], end = xtspar[2], deltat = xtspar[3])
}
"traceplot" <-
function (x, combine.chains = F, smooth = T, ...) 
{
        x <- as.mcmc(x)
        if (combine.chains) 
                for (i in 1:nvar(x)) matplot(time(x), x[, i, 
                        , drop = T], type = "l", xlab = "Iterations", 
                        ylab = varnames(x)[i])
        else for (i in 1:nchain(x)) {
                for (j in 1:nvar(x)) {
                        xp <- as.vector(time(x))
                        yp <- x[, j, i, drop = T]
                        plot(xp, yp, type = "l", xlab = "Iterations", 
                                ylab = varnames(x)[j])
                        if (smooth) 
                                lines(lowess(xp, yp), col = "red")
                        mtext(paste("Trace of", varnames(x)[j]), 
                                side = 3, line = 2.5)
                        mtext(paste("(", format(niter(x)), " values per trace)", 
                                sep = ""), side = 3, line = 0.65)
                }
        }
}
"tspar" <-
function (x) 
attr(x, "tspar")
"tspar<-" <-
function (x, value) 
{
        attr(x, "tspar") <- value
        x
}
"window.mcmc" <-
function (x, start, end, thin) 
{
        x <- as.mcmc(x)
        xtspar <- tspar(x)
        xstart <- xtspar[1]
        xend <- xtspar[2]
        xthin <- xtspar[3]
        if (missing(thin)) 
                thin <- xthin
        else if (thin%%xthin != 0) {
                thin <- xthin
                warning("Thin value not changed")
        }
        xtime <- as.vector(time(x))
        ts.eps <- .Options$ts.eps
        if (missing(start)) 
                start <- xstart
        else if (length(start) != 1) 
                stop("bad value for start")
        else if (start < xstart) {
                start <- xstart
                warning("start value not changed")
        }
        if (missing(end)) 
                end <- xend
        else if (length(end) != 1) 
                stop("bad value for end")
        else if (end > xend) {
                end <- xend
                warning("end value not changed")
        }
        if (start > end) 
                stop("start cannot be after end")
        if (all(abs(xtime - start) > abs(start) * ts.eps)) {
                start <- xtime[(xtime > start) & ((start + xthin) > 
                        xtime)]
        }
        if (all(abs(end - xtime) > abs(end) * ts.eps)) {
                end <- xtime[(xtime < end) & ((end - xthin) < 
                        xtime)]
        }
        use <- 1:nrow(x)
        use <- use[use >= trunc((start - xstart)/xthin + 1.5) & 
                use <= trunc((end - xstart)/xthin + 1.5) & (use - 
                trunc((start - xstart)/xthin + 1.5))%%(thin%/%xthin) == 
                0]
        return(x[use, , ])
}
