CODA _ function() { # # CODA v0.40bs1 -- 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) # Brian Smith (University of Iowa, Iowa City, US) # # Date of last update: 2/13/99 # Edited by BS to fix incompatibilities with # S-PLUS v5.0 # Added the "globals.s" library to centralized the # the management of all globals variables # options("object.size"=1e+9) globals.create() globals.init() on.exit(tidy.up()) home<-"" if(is.null(globals()$parmsmat)) { cat("\n _______________________________________________________________ \n|\t\t\t\t\t\t\t\t|\n|\t\t Welcome to CODA!\t\t\t\t|\n| Convergence Diagnostics and Output Analysis for BUGS output |\n|_______________________________________________________________|\n|\t\t\t\t\t\t\t\t|\n|\t\t\t\t\t\t\t\t|\n") cat("| Authors : Nicky Best, Kate Cowles & Karen Vines. \t\t|\n|\t\t\t\t\t\t\t\t|\n|\t\t\t\t\t\t\t\t|\n") cat("| CODA : Copyright (c) 1995 MRC Biostatistics Unit.\t\t|\n| All rights reserved\t\t\t\t\t\t|\n| Version 0.40bs1\t\t\t\t\t\t\t|\n|\t\t\t\t\t\t\t\t|\n") cat("|_______________________________________________________________|\n") # Ask user whether they are running S-Plus under UNIX or WINDOWS cat("\nAre you are running S-Plus under WINDOWS (y or n)?:\n") win <- scan(what = character(), n=1, sep = "\n", strip.white=T) if(win=="y" | win=="Y") { home <- set.windows.home() } globals(Home = home) globals(outfile = paste(globals()$Home, "CODA.LOG", sep="")) cat("\n _______________________________________________________________ \n|\t\t\t\t\t\t\t\t|\n|\t\t Welcome to CODA!\t\t\t\t|\n| Convergence Diagnostics and Output Analysis for BUGS output |\n|_______________________________________________________________|\n|\t\t\t\t\t\t\t\t|\n|\t\t\t\t\t\t\t\t|\n", file=globals()$outfile, append=F) cat("| Authors : Nicky Best, Kate Cowles & Karen Vines. \t\t|\n|\t\t\t\t\t\t\t\t|\n|\t\t\t\t\t\t\t\t|\n", file=globals()$outfile, append=F) cat("| CODA : Copyright (c) 1995 MRC Biostatistics Unit.\t\t|\n| All rights reserved\t\t\t\t\t\t|\n| Version 0.30\t\t\t\t\t\t\t|\n|\t\t\t\t\t\t\t\t|\n", file=globals()$outfile, append=T) cat("|_______________________________________________________________|\n\n", file=globals()$outfile, append=T) cat("\n\nDo you wish to:\n") file.menu <- c( "Begin a new CODA session using BUGS output files", "Begin a new CODA session using data saved from a previous session", "Quit") pick <- menu(file.menu)+1 } else { cat("\n\nDo you wish to:\n") file.menu <- c("Continue with current CODA session", "Begin a new CODA session using BUGS output files", "Begin a new CODA session using data saved from a previous session", "Quit") pick <- menu(file.menu) } switch(pick, main(), { cat("\n\nEnter BUGS output filenames, separated by return key\n(Leave blank and hit return to exit):\n") globals(filenames = scan(what = character(), sep = "\n", strip.white=T)) }, { cat("\nEnter name of S-Plus data object saved from a previous CODA session:\n") while(T) { outname <- scan(what = character(), sep = "\n", n=1, strip.white=T) if(outname==0) { break } if(exists(outname)) { local.temp <- eval(parse(text = paste("", outname, sep = ""))) is.CODA <- F if(is.list(local.temp)) { if(is.null(names(local.temp))) { names(local.temp) <- paste("chain_", 1:length(local.temp), sep="") } is.mtrx <- numeric(length(local.temp)) dim.mtrx <- matrix(nrow=length(local.temp), ncol=2) for(i in 1:length(local.temp)) { is.mtrx[i] <- is.matrix(local.temp[[i]]) if(is.mtrx[i]) dim.mtrx[i,] <- dim(local.temp[[i]]) } if(all(is.mtrx)) { if(all(dim.mtrx[,1]==dim.mtrx[1,1]) & all(dim.mtrx[,2]==dim.mtrx[1,2])) { is.CODA <- T } } globals(out = local.temp) } if(is.CODA) { break } else { cat("\n", outname, " is not a CODA object\n", sep="") } } else { cat("\nCan't find this file!\n") } cat("Please re-enter another name or type 0 to quit:\n") } if(outname==0) { globals(out = NULL) return(eval(parse(text="cat('')"))) break } globals(filenames = names(globals()$out)) globals(Labels = dimnames(globals()$out[[1]])[[2]][-1]) globals(nfiles = length(globals()$filenames)) }, { return(eval(parse(text="cat('')"))) break } ) if(globals()$quit.CODA) { break } cat("\nEnter problem title, followed by return key:\n") globals(mtitle = scan(what = character(), sep = "\n", n = 1, strip.white=T)) cat("Problem title:", globals()$mtitle, "\n\n", file=globals()$outfile, append=T) if(pick == 2) { globals(nfiles = length(globals()$filenames)) indices <- list() outtmp <- list() niters <- numeric() start.iter <- numeric() restart(TRUE) for(i in 1:globals()$nfiles) { indices[[i]] <- inddat(file = paste(globals()$filenames[i], ".ind", sep = "")) outtmp[[i]] <- readdat(vnames = indices[[i]][[1]], index = indices[[i]], file = paste(globals()$filenames[i], ".out", sep = "")) niters[i] <- nrow(outtmp[[i]]) start.iter[i] <- outtmp[[i]][1, "iter"] } # # If > 1 input file, check that each file is of same length, starts # from same iteration, and contains the same number of variables # if(globals()$nfiles > 1) { parm.OK <- rep(T, ncol(outtmp[[1]])-1) for(i in 1:(ncol(outtmp[[1]])-1)) { j <- 2 while(parm.OK[i] && j <= globals()$nfiles) { parm.OK[i]<-ifelse(any(indices[[j]]$varname == indices[[1]]$varname[i]), T, F) j <- j+1 } } if(all(parm.OK == F)) { cat("\n******* Warning: *******\nThere are no variables common to all chains\nOnly first chain (", globals()$filenames[1], ") will be used\n", sep="") length(outtmp) <- 1 length(indices) <- 1 local.temp <- globals()$filenames length(local.temp) <- 1 globals(filenames = local.temp) globals(nfiles = 1) } else { if(any(parm.OK == F)) { cat("\n******* Warning: *******\nOnly variable(s) ") cat(indices[[1]]$varname[parm.OK], sep=", ") cat(" are common to all chains\nAll others will be deleted\n") for(i in 1:globals()$nfiles) { outtmp[[i]] <- outtmp[[i]][, c("iter", indices[[1]]$varname[parm.OK])] } tmpind <- indices[[1]]$varname[parm.OK] for(i in 1:globals()$nfiles) { indices[[i]]$varname <- tmpind } } if(!all(niters == niters[1])) { miniter <- min(niters) cat("\n******* Warning: *******\nChains are of different lengths\nTruncating each to last", miniter, "iterations\n") for(i in 1:globals()$nfiles) { outtmp[[i]] <- outtmp[[i]][(niters[i] - miniter+1):niters[i], ] } } if(!all(start.iter == start.iter[1])) { cat("\n******* Warning: *******\nChains start at different iteration numbers\nAll plots will be labelled by iteration numbers for first chain (", globals()$filenames[1], ")\n", sep="") } } } # # Now sort columns of "outtmp" so that variables for each chain # appear in the same order # local.temp <- list() for(i in 1:globals()$nfiles) { local.temp[[i]] <- matrix(ncol = ncol(outtmp[[i]]), nrow = nrow(outtmp[[i]])) local.temp[[i]][, 1] <- outtmp[[i]][, 1] local.temp[[i]][, -1] <- outtmp[[i]][, order(indices[[i]]$varname) + 1] indices[[i]]$varname <- sort(indices[[i]]$varname) dimnames(local.temp[[i]]) <- list(NULL, c("iter", indices[[i]]$varname)) } names(local.temp) <- globals()$filenames globals(out = local.temp) globals(Labels = indices[[1]]$varname) } outtmp <- globals()$out globals(nparms = ncol(globals()$out[[1]]) - 1) # # Format labels and filenames for printing purposes: # globals(labels.list = pretty.print(matrix(c("VARIABLE NUMBER", 1:length(globals()$Labels), "VARIABLE NAME", globals()$Labels), ncol=2), print.it=F)) globals(short.labels = abbreviate(globals()$Labels, minlength=7)) globals(Labels = abbreviate(globals()$Labels, minlength=12)) globals(files.list = pretty.print(matrix(c("CHAIN NUMBER", 1:globals()$nfiles, "CHAIN NAME", globals()$filenames), ncol=2), print.it=F)) # # Set the data.defaults which depend on the working data # data.defaults(parms = 1:globals()$nparms, chains = 1:globals()$nfiles, end = dim(globals()$out[[1]])[1]) # # Prompt user to specify which (if any) variables take values # > 0 or lie in the interval (0, 1). These will be transformed # to the log or logit scale respectively for Gelman & Rubin's # diagnostic; for kernel density estiamtes, the negative reflection # technique will be used to ensire estiamtes within the appropriate # range # logit.var() log.var() # # Create working data array containing only those samples # selected under the current data default settings # globals(first.time = T) cat("\nCreating working data object...\n") working.data() globals(first.time = F) # # Call main menu # main() } # # ----------------------------------------------------------------------- # set.windows.home _ function() { # # set.windows.home -- allows WINDOWS users to change home directory # # Author: Nicky Best # Date of last update: 28/8/97 # # Ask user to define working directory tmp <- search()[1] wkdir <- substring(tmp, 1, nchar(tmp)-6) slash<-0 for(i in 1:nchar(wkdir)) { if(substring(wkdir,i,i) == "\\") slash<-append(slash, i) } wkdir2<-NULL for(j in 1: (length(slash)-1)) { wkdir2 <- paste(wkdir2, paste(substring(wkdir, slash[j]+1, slash[j+1]), "\\", sep=""), sep="") } wkdir2 <- paste(wkdir2, substring(wkdir,slash[length(slash)]+1, nchar(wkdir)), sep="") cat("\n\nThe current working directory is", wkdir2, "\n\nBy default, CODA will assume all input files (i.e. *.out \nand *.ind files) are located in this directory,\nand will write the CODA.LOG file to this directory\n\nTo change the working directory, enter the new \ndirectory name and full path, or leave blank and \nhit return to continue.\n(Remember to use a DOUBLE backslash", paste("\\","\\", sep=""), "to separate \ndirectory names e.g.", paste("C:", "\\", "\\", "SPLUSWIN)", sep=""), "\n") home <- scan(what = character(), n=1, sep = "\n", strip.white=T) if(length(home)>0) { home <- paste(home, "\\", sep= "") } else { home <- paste(wkdir, "\\", sep= "") } home } # # ----------------------------------------------------------------------- # logit.var _ function() { # # logit.var -- Prompts user to specify variables which lie between (0, 1) # # Author: Nicky Best # while(T) { cat("\nAre any variables restricted to values between 0 and 1 (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } if(ans == "Y" | ans == "y") { cat("\nAvailable variables:\n\n") cat(globals()$labels.list) while(T) { cat("\nEnter relevant variable number(s), separated by commas \n(Ranges such as 3:7 may be specified)\n(Enter 0 for none):\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { func.defaults(prop.parms = numeric()) for(i in 1:length(ans)) { func.defaults(prop.parms = append(func.defaults()$prop.parms, eval(parse(text = ans[i])))) } if(min(func.defaults()$pos.parms) < 0 | max(func.defaults()$prop.parms) > length(globals()$Labels)) { cat("Error: you have specified variable number(s) out of range!\n\n") } else break } } } # # Check that pos.parms and prop.parms do not include same variables # repeat{ tmp <- rep(T, length(func.defaults()$pos.parms)) if(length(func.defaults()$pos.parms) > 0 && length(func.defaults()$prop.parms) > 0) { for(i in 1:length(func.defaults()$pos.parms)) { tmp[i] <- ifelse(any(func.defaults()$prop.parms == func.defaults()$pos.parms[i]), T, F) } } if(any(tmp) && length(func.defaults()$prop.parms) > 0 && func.defaults()$prop.parms[1] != 0) { cat("\nError: You have specified variable(s) ") cat(globals()$Labels[func.defaults()$pos.parms[tmp]], sep=", ") cat(" to be BOTH positive and on (0, 1)\nOnly ONE of these options is allowed per variable\n") while(T) { cat("Please re-enter number(s) of variables taking positive values:\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { func.defaults(pos.parms = numeric()) for(i in 1:length(ans)) { func.defaults(pos.parms = append(func.defaults()$pos.parms, eval(parse(text = ans[i])))) } if(min(func.defaults()$pos.parms) < 0 | max(func.defaults()$pos.parms) > length(globals()$Labels)) { cat("Error: you have specified variable number(s) out of range!\n\n") } else break } } } while(T) { cat("Now re-enter number(s) of variables taking values on (0, 1):\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { func.defaults(prop.parms = numeric()) for(i in 1:length(ans)) { func.defaults(prop.parms = append(func.defaults()$prop.parms, eval(parse(text = ans[i])))) } if(min(func.defaults()$prop.parms) < 0 | max(func.defaults()$prop.parms) > length(globals()$Labels)) { cat("Error: you have specified variable number(s) out of range!\n\n") } else break } } } } else break } changes(pos.parms = T) } } # # ----------------------------------------------------------------------- # log.var _ function() { # # log.var -- Prompts user to specify all-positive variables # # Author: Nicky Best # while(T) { cat("\nAre any variables restricted to all positive values (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } if(ans == "Y" | ans == "y") { cat("\nAvailable variables:\n\n") cat(globals()$labels.list) while(T) { cat("\nEnter relevant variable number(s), separated by commas \n(Ranges such as 3:7 may be specified)\n(Enter 0 for none):\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { func.defaults(pos.parms = numeric()) for(i in 1:length(ans)) { func.defaults(pos.parms = append(func.defaults()$pos.parms, eval(parse(text = ans[i])))) } if(min(func.defaults()$pos.parms) < 0 | max(func.defaults()$pos.parms) > length(globals()$Labels)) { cat("Error: you have specified variable number(s) which are out of range!\n") } else break } } } # # Check that pos.parms and prop.parms do not include same variables # repeat{ tmp <- rep(T, length(func.defaults()$pos.parms)) if(length(func.defaults()$pos.parms) > 0 && length(func.defaults()$prop.parms) > 0) { for(i in 1:length(func.defaults()$pos.parms)) { tmp[i] <- ifelse(any(func.defaults()$prop.parms == func.defaults()$pos.parms[i]), T, F) } } if(any(tmp) && length(func.defaults()$prop.parms) > 0 && func.defaults()$prop.parms[1] != 0) { cat("\nError: You have specified variable(s) ") cat(globals()$Labels[func.defaults()$pos.parms[tmp]], sep=", ") cat(" to be BOTH positive and on (0, 1)\nOnly ONE of these options is allowed per variable\n") while(T) { cat("Please re-enter number(s) of variables taking positive values:\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { func.defaults(pos.parms = numeric()) for(i in 1:length(ans)) { func.defaults(pos.parms = append(func.defaults()$pos.parms, eval(parse(text = ans[i])))) } if(min(func.defaults()$pos.parms) < 0 | max(func.defaults()$pos.parms) > length(globals()$Labels)) { cat("Error: you have specified variable number(s) out of range!\n\n") } else break } } } while(T) { cat("Now re-enter number(s) of variables taking values on (0, 1):\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { func.defaults(prop.parms = numeric()) for(i in 1:length(ans)) { func.defaults(prop.parms = append(func.defaults()$prop.parms, eval(parse(text = ans[i])))) } if(min(func.defaults()$prop.parms) < 0 | max(func.defaults()$prop.parms) > length(globals()$Labels)) { cat("Error: you have specified variable number(s) out of range!\n\n") } else break } } } } else break } changes(pos.parms = T) } } # # ----------------------------------------------------------------------- # main _ function() { # # main -- CODA main menu # # Author: Kate Cowles # Modified by: Nicky Best # # Note that using the format: # # while(T) { # this.menu <- .... # pick <- menu(this.menu) # switch(pick, # ..., # break # ) # } # # doesn't work here, as main() may have been called from another # menu e.g. change.defaults, in which case, break will just return # control to that menu, rather than quitting. Hence the need for # the logical `quit.CODA', and the `while(pick != x)' in the # other menus, so that they do not get re-evaluated after the # `Quit' option has been selected from the main coda.menu # restart(T) globals(pick = 0) globals(quit.CODA = F) while(globals()$pick != 4) { if(globals()$quit.CODA) break cat("\nCODA Main Menu\n**************\n") coda.menu <- c("Output Analysis", "Diagnostics", "List/Change Defaults", "Quit") globals(pick = menu(coda.menu)) switch(globals()$pick, outanal(), diags(), change.defaults(), { globals(quit.CODA = T) break } ) } } # # ----------------------------------------------------------------------- # diags _ function() { # # diags -- Diagnostics menu for CODA system # # Author: Kate Cowles # restart(T) globals(pick1 = 0) while(globals()$pick1 != 7 && globals()$pick1 != 8) { cat("\nCODA Diagnostics Menu\n*********************\n") diag.menu <- c("Geweke", "Gelman and Rubin", "Raftery and Lewis", "Heidelberger and Welch", "Autocorrelations", "Cross-Correlations", "List/Change Defaults", "Return to Main Menu") globals(pick1 = menu(diag.menu)) switch(globals()$pick1, geweke(), gandr(), randl2(), heidel(), autocorr(), parmcorr(), change.defaults(), main() ) } } # # ----------------------------------------------------------------------- # outanal _ function() { # # outanal -- Output analysis menu for CODA system # # Author: Kate Cowles # restart(T) globals(pick2 = 0) while(globals()$pick2 != 3 && globals()$pick2 != 4) { cat("\nCODA Output Analysis Menu\n*************************\n") diag.menu <- c("Plots", "Statistics", "List/Change Defaults", "Return to Main Menu") globals(pick2 = menu(diag.menu)) switch(globals()$pick2, outplots(), outstats(), change.defaults(), main() ) } } # # ----------------------------------------------------------------------- # change.defaults _ function() { # # change.defaults -- Defaults menu for CODA system # # Author: Nicky Best # restart(T) globals(pick3 = 0) while(globals()$pick3 != 9 && globals()$pick3 != 10 && globals()$pick3 != 11) { cat("\nCODA Main Defaults Menu\n***********************\n") defaults.menu <- c("List current defaults", "Select variables for analysis", "Select chains for analysis", "Select iterations for analysis", "Select thinning interval", "Plots - Defaults Menu", "Summary Statistics - Defaults Menu", "Diagnostics - Defaults Menu", "Return to Output Analysis Menu", "Return to Diagnostics Menu", "Return to Main Menu") globals(pick3 = menu(defaults.menu)) switch(globals()$pick3, { # # List current defaults # display.defaults(data=T, funcs=T) }, { # # Select variables for analysis # cat("\nAvailable variables:\n\n") cat(globals()$labels.list) while(T) { cat("\nEnter variable number(s) you wish to analyse, separated by commas \n(Ranges such as 3:7 may be specified):\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { data.defaults(parms = numeric()) for(i in 1:length(ans)) { data.defaults(parms = append(data.defaults()$parms, eval(parse(text = ans[i])))) } if(min(data.defaults()$parms) < 1 | max(data.defaults()$parms) > length(globals()$Labels)) { cat("Error: you have specified variable number(s) which are out of range!\n") } else { data.defaults(parmsC = numeric()) for(i in 1:globals()$nparms) { if(all(data.defaults()$parms != i)) data.defaults(parmsC = append(data.defaults()$parmsC, i)) } break } } } } changes(parms = T) }, { # # Select chains for analysis # cat("\nAvailable chains:\n\n") cat(globals()$files.list) while(T) { cat("\nEnter chain number(s) you wish to analyse, separated by commas \n(Ranges such as 3:7 may be specified):\n") ans <- scan(what = character(), sep = ",", strip.white=T) if(length(ans) > 0) { err <- logical() for(i in 1:length(ans)) { err[i] <- any(is.na(as.numeric((substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))[ substring(ans[i], 1:nchar(ans[i]), 1:nchar(ans[i]))!=":"])))) } if(any(err)) { cat("Error: you have specified a non-numeric value!\n") } else { data.defaults(chains = numeric()) for(i in 1:length(ans)) { data.defaults(chains = append(data.defaults()$chains, eval(parse(text = ans[i])))) } if(min(data.defaults()$chains) < 1 | max(data.defaults()$chains) > length(globals()$filenames)) { cat("Error: you have specified chain number(s) which are out of range!\n") } else { data.defaults(chainsC = numeric()) for(i in 1:globals()$nfiles) { if(all(data.defaults()$chains != i)) data.defaults(chainsC = append(data.defaults()$chainsC, i)) } break } } } } changes(chains = T) }, { # # Select iterations for analysis # while(T) { cat("\nIterations available = ", globals()$out[[1]][1, "iter"], ":", dim(globals()$out[[1]])[1]+globals()$out[[1]][1, "iter"]-1, "\n", "\n\nEnter iteration you wish to start at:\n", sep="") data.defaults(start = scan(what = numeric(), n = 1)) if(length(data.defaults()$start) > 0) break } while(T) { cat("Enter iteration you wish to end at:\n") data.defaults(end = scan(what = numeric(), n = 1)) if(length(data.defaults()$end) > 0) break } # # Check that end >= start and iterations in range # while(data.defaults()$end < data.defaults()$start) { while(T) { cat("\nError: End is less than start\nPlease re-enter starting iteration:\n") data.defaults(start = scan(what = numeric(), n = 1)) if(length(data.defaults()$start) > 0) break } while(T) { cat("Now re-enter end iteration:\n") data.defaults(end = scan(what = numeric(), n = 1)) if(length(data.defaults()$end) > 0) break } } while(data.defaults()$start < globals()$out[[1]][1, "iter"] | data.defaults()$end > (dim(globals()$out[[1]])[1] +globals()$out[[1]][1, "iter"])) { while(T) { cat("\nError: Iterations requested are out of range\nPlease re-enter starting iteration:\n") data.defaults(start = scan(what = numeric(), n = 1)) if(length(data.defaults()$start) > 0) break } while(T) { cat("Now re-enter end iteration:\n") data.defaults(end = scan(what = numeric(), n = 1)) if(length(data.defaults()$end) > 0) break } } data.defaults(start = data.defaults()$start-globals()$out[[1]][1, "iter"]+1) data.defaults(end = data.defaults()$end-globals()$out[[1]][1, "iter"]+1) while(T) { cat("Do you wish to retain ALL iterations for trace plots:\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } data.defaults(all.trace = ifelse((ans == "Y" | ans == "y"), T, F)) changes(iter = T) }, { # # Thinning interval # while(T) { cat("\nEnter thinning interval:\n") data.defaults(thin = scan(what = numeric(), n = 1)) if(length(data.defaults()$thin) > 0) break } changes(thin = T) }, plot.defaults(), stats.defaults(), diag.defaults(), { if(any(c(changes()$parms, changes()$chains, changes()$iter, changes()$thin, changes()$pos.parms))) { # # Re-create working data array: # working.data() } outanal() }, { if(any(c(changes()$parms, changes()$chains, changes()$iter, changes()$thin, changes()$pos.parms))) { # # Re-create working data array: # working.data() } diags() }, { if(any(c(changes()$parms, changes()$chains, changes()$iter, changes()$thin, changes()$pos.parms))) { # # Re-create working data array: # working.data() } main() } ) } } # # ----------------------------------------------------------------------- # plot.defaults _ function() { # # plot.defaults -- Plotting defaults menu for CODA system # # Author: Nicky Best # restart(T) while(T) { cat("\nCODA Plotting Defaults Menu\n***************************\n") plot.defaults.menu <- c("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", "Exponent for scientific notation", "Return to Previous Menu") pick4 <- menu(plot.defaults.menu) switch(pick4, { # # Plot trace of samples # while(T) { cat("\nDo you want plots of variable traces (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(trace = ifelse((ans == "Y" | ans == "y"), T, F)) }, { # # Plot kernel density estimate # while(T) { cat("\nDo you want density plots (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(densplot = ifelse((ans == "Y" | ans == "y"), T, F)) }, { # # Add smooth line through trace plot # while(T) { cat("\nDo you want a smoothed line superimposed on the trace plots (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(lowess = ifelse((ans == "Y" | ans == "y"), T, F)) }, { # # Select separate plots per chain # while(T) { cat("\nDo you want separate plots per chain for each variable (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(sepplot = ifelse((ans == "Y" | ans == "y"), T, F)) }, { # # Select separate page of plots per variable # while(T) { cat("\nDo you want plots for each variable to appear on separate pages (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(onepage = ifelse((ans == "Y" | ans == "y"), T, F)) }, { # # Specify page layout for plots # while(T) { cat("\nDo you want to specify your own page layout for the plots (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(user.layout = ifelse((ans == "Y" | ans == "y"), T, F)) if(func.defaults()$user.layout) { while(T) { cat("Enter number of rows of plots per page (maximum=7):\n") func.defaults(mrows = scan(what = numeric(), n = 1)) if(length(func.defaults()$mrows) > 0) break } while(T) { cat("Enter number of columns of plots per page (maximum=8):\n") func.defaults(mcols = scan(what = numeric(), n = 1)) if(length(func.defaults()$mcols) > 0) break } # # Check that mrows x mcols <= 7 x 8 # while(func.defaults()$mrows > 7 | func.defaults()$mcols > 8) { while(T) { cat("\nError: maximum of 7 rows and 8 columns allowed\nPlease re-enter number of rows:\n") func.defaults(mrows = scan(what = numeric(), n = 1)) if(length(func.defaults()$mrows) > 0) break } while(T) { cat("\nNow re-enter number of columns:\n") func.defaults(mcols = scan(what = numeric(), n = 1)) if(length(func.defaults()$mcols) > 0) break } } } }, kernel.width(), ps.file.options(), { # # Scientific notation # while(T) { cat("\nEnter exponent for scientific notation\n(Default = 6 i.e. numbers with exponents < -6 or > 6 will be printed in scientific notation):\n") func.defaults(scientific = round(scan(what = numeric(), n = 1))) if(length(func.defaults()$scientific) > 0) break } }, break) } } # # ----------------------------------------------------------------------- # kernel.width _ function() { # # kernel.width -- Function called by plotting defaults menu to allow user # to change bandwidth function used for kernel density estimation. # This function is called separately to allow inclusion of the # `restart(T)' function - this is an S-Plus function which # causes a function to be restarted if errors occur, rather than # dumping and return control to the top level. If the # user selects option 3 (i.e. user defined bandwidth) it is easy # to miss-type the function they wish to specify, leading to syntax # errors when the function is parsed, which would otherwise cause # CODA to crash # # Author: Nicky Best # restart(T) # # Select bandwidth function for kernel smoothing # if(!func.defaults()$densplot) { cat("\nNo density plots requested - this option is irrelevant\n") } else { cat("\nSelect kernel bandwidth function required:\n") kernel.menu <- c("Default: Smooth (0.25 * sample range)", "Coarse (Silverman 1986 eqn. 3.28 & 3.30)", "User-defined function", "Return to Plotting Defaults Menu") pick1 <- menu(kernel.menu) switch(pick1, func.defaults(bandwidth = function(y) { (max(y)-min(y))/4 }), func.defaults(bandwidth = function(y) { 1.06 * min(sqrt(var(y)), (quantile(y,0.75)-quantile(y,0.25)) /1.34)*(length(y)^(-1/5)) }), { # # Allow user to enter thir 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 # 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) { func.defaults(bandwidth = "function(y){") for(i in 1:length(ans)) { func.defaults(bandwidth = paste(func.defaults()$bandwidth, ans[i], sep="")) } func.defaults(bandwidth = paste( func.defaults()$bandwidth, "}", sep="")) func.defaults(bandwidth = eval(parse(text = func.defaults()$bandwidth))) # # Carry out simple test to check whether the # function entered makes sense # func.OK <- test.bandwidth() } } }, plot.defaults() ) } } # # ----------------------------------------------------------------------- # 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 <- func.defaults()$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 } # # ----------------------------------------------------------------------- # ps.file.options _ function() { # # ps.file.options -- Function called by plotting defaults menu to allow user # to select size and orientation of any plots to be saved # as a postscript file # # Author: Nicky Best # restart(T) cat("\nSelect options required when saving plots to postscript files:\n") ps.menu <- c("Default: Full page portrait", "Full page landscape", "Half page portrait") func.defaults(ps.plot = menu(ps.menu)) } # # ----------------------------------------------------------------------- # stats.defaults _ function() { # # stats.defaults -- Summary statistics defaults menu for CODA system # # Author: Nicky Best # restart(T) while(T) { cat("\nCODA Summary Statistics Defaults Menu \n*************************************\n") stats.defaults.menu <- c("Combine chains for summary statistics", "Batch size for calculating batch means", "Quantiles for summary statistics", "Number of significant digits for printing", "Return to Previous Menu") pick4 <- menu(stats.defaults.menu) switch(pick4, { # # Pool over chains for summary stats # while(T) { cat("\nDo you want combine all chains to calculate summary statistics (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(combine.stats = ifelse((ans == "Y" | ans == "y"), T, F)) }, { # # Batch size for calculating batch means # while(T) { cat("\nEnter batch size (Default = 25):\n") func.defaults(batch.size = scan(what = numeric(), n = 1)) if(length(func.defaults()$batch.size) > 0) break } }, { # # Quantiles # while(T) { cat("\nEnter quantiles required, separated by commas\n(Default = 0.025, 0.25, 0.5, 0.75, 0.975):\n") func.defaults(quantiles = scan(what = numeric(), sep=",")) if(length(func.defaults()$quantiles) > 0) break } # # Check that quantiles are between 0 and 1 # while(any(func.defaults()$quantiles < 0) || any(func.defaults()$quantiles > 1)) { while(T) { cat("\nError: Quantiles out of range (0, 1)\nPlease re-enter required quantiles:\n") func.defaults(quantiles = scan(what = numeric(), sep=",")) if(length(func.defaults()$quantiles) > 0) break } } }, { # # Significant digits # while(T) { cat("\nEnter number of significant digits to be printed (Default = 3):\n") func.defaults(digits = round(scan(what = numeric(), n = 1))) if(length(func.defaults()$digits) > 0) break } }, break ) } } # # ----------------------------------------------------------------------- # diag.defaults _ function() { # # diag.defaults -- Diagnostics defaults menu for CODA system # # Author: Nicky Best # restart(T) while(T) { cat("\nCODA Diagnostics Defaults Menu\n******************************\n") diag.defaults.menu <- c("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 Previous Menu") pick4 <- menu(diag.defaults.menu) switch(pick4, { # # Select window widths for Geweke's diagnostic # while(T) { cat("\nEnter fraction of chain to include in 1st window of \nGeweke's diagnostic (Default = 0.1):\n") func.defaults(frac1 = scan(what = numeric(), n = 1)) if(length(func.defaults()$frac1) > 0) break } while(T) { cat("\nEnter fraction of chain to include in 2nd window of \nGeweke's diagnostic (Default = 0.5):\n") func.defaults(frac2 = scan(what = numeric(), n = 1)) if(length(func.defaults()$frac2) > 0) break } # # Check that sum of fractions doesn't exceed 1.0 # while(func.defaults()$frac1+func.defaults()$frac2 >= 1) { while(T) { if(func.defaults()$frac1+func.defaults()$frac2 == 1) { cat("\nError: 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:\n") } else { cat("\nError: Sum of fractions in 1st and 2nd windows exceeds 1.0\nPlease re-enter fraction of chain to include in 1st window:\n") } func.defaults(frac1 = scan(what = numeric(), n = 1)) if(length(func.defaults()$frac1) > 0) break } while(T) { cat("\nNow re-enter fraction of chain to include in 2nd window:\n") func.defaults(frac2 = scan(what = numeric(), n = 1)) if(length(func.defaults()$frac2) > 0) break } } }, { # # Bin size for plotting Geweke's daignostic # cat("\nOptions for defining bin size to plot Geweke's diagnostic:\n") geweke.menu <- c("Default: bin width = 10; maximum number of bins = 50", "User-specified bin width", "User-specified total number of bins") pick1 <- menu(geweke.menu) switch(pick1, { func.defaults(geweke.max = 50) func.defaults(geweke.bin = 10) }, { func.defaults(geweke.max = F) while(T) { cat("\nEnter required bin width (Default = 10 iterations):\n") func.defaults(geweke.bin = scan(what = numeric(), n = 1)) if(length(func.defaults()$geweke.bin) > 0) break } # # Check that bin width is between 1 and Niter-50 # while(func.defaults()$geweke.bin < 1 || func.defaults()$geweke.bin > (globals()$Niters-50)) { while(T) { cat("\nError: Bin width out of range\nPlease re-enter a value between 1 and ", globals()$Niters-50, ":\n", sep="") func.defaults(geweke.bin=scan(what = numeric(), n = 1)) if(length(func.defaults()$geweke.bin) > 0) break } } }, { func.defaults(geweke.bin = F) while(T) { cat("\nEnter total number of bins required (Default = 50):\n") func.defaults(geweke.max = scan(what = numeric(), n = 1)) if(length(func.defaults()$geweke.max) > 0) break } # # Check that maximum is between 2 and Niter-50 # while(func.defaults()$geweke.max < 2 || func.defaults()$geweke.max > (globals()$Niters-50)) { while(T) { cat("\nError: Number out of range\nPlease re-enter a value between 2 and ", globals()$Niters-50, ":\n", sep="") func.defaults(geweke.max=scan(what = numeric(), n = 1)) if(length(func.defaults()$geweke.max) > 0) break } } func.defaults(geweke.max = func.defaults()$geweke.max - 1) # (need to subtract 1 as first bin of size 50 is automatic) } ) }, { # # Bin size for plotting Gelman & Rubin's daignostic # cat("\nOptions for defining bin size to plot Gelan & Rubin's diagnostic:\n") gr.menu <- c("Default: bin width = 10; maximum number of bins = 50", "User-specified bin width", "User-specified total number of bins") pick1 <- menu(gr.menu) switch(pick1, { func.defaults(gr.max = 50) func.defaults(gr.bin = 10) }, { func.defaults(gr.max = F) while(T) { cat("\nEnter required bin width (Default = 10 iterations):\n") func.defaults(gr.bin = scan(what = numeric(), n = 1)) if(length(func.defaults()$gr.bin) > 0) break } # # Check that bin width is between 1 and Niter-50 # while(func.defaults()$gr.bin < 1 || func.defaults()$gr.bin > (globals()$Niters-50)) { while(T) { cat("\nError: Bin width out of range\nPlease re-enter a value between 1 and ", globals()$Niters-50, ":\n", sep="") func.defaults(gr.bin=scan(what = numeric(), n = 1)) if(length(func.defaults()$gr.bin) > 0) break } } }, { func.defaults(gr.bin = F) while(T) { cat("\nEnter total number of bins required (Default = 50):\n") func.defaults(gr.max = scan(what = numeric(), n = 1)) if(length(func.defaults()$gr.max) > 0) break } # # Check that maximum is between 2 and Niter-50 # while(func.defaults()$gr.max < 2 || func.defaults()$gr.max > (globals()$Niters-50)) { while(T) { cat("\nError: Number out of range\nPlease re-enter a value between 2 and ", globals()$Niters-50, ":\n", sep="") func.defaults(gr.max=scan(what = numeric(), n = 1)) if(length(func.defaults()$gr.max) > 0) break } } func.defaults(gr.max = func.defaults()$gr.max - 1) # (need to subtract 1 as first bin of size 50 is automatic) } ) }, { # # Select parameters for Raftery & Lewis' diagnostic # while(T) { cat("\nEnter quantile to be estimated for Raftery & Lewis' diagnostic\n(Default = 0.025):\n") func.defaults(q = scan(what = numeric(), n = 1)) if(length(func.defaults()$q) > 0) break } # # Check that quantile is between 0 & 1 # while(func.defaults()$q < 0 || func.defaults()$q > 1) { while(T) { cat("\nError: Quantile out of range\nPlease re-enter a value between 0 and 1:\n") func.defaults(q = scan(what = numeric(), n = 1)) if(length(func.defaults()$q) > 0) break } } while(T) { cat("\nEnter required precision (Default = 0.005):\n") func.defaults(r = scan(what = numeric(), n = 1)) if(length(func.defaults()$r) > 0) break } # # Check that precision is sensible # if(func.defaults()$r > func.defaults()$q) { while(T) { cat("\nWarning: You have specified a very wide precision interval\nPlease re-enter this value to confirm, or a different value to change:\n") func.defaults(r = scan(what = numeric(), n = 1)) if(length(func.defaults()$r) > 0) break } } while(T) { cat("\nEnter required probability (Default = 0.95):\n") func.defaults(s = scan(what = numeric(), n = 1)) if(length(func.defaults()$s) > 0) break } # # Check that probability is between 0 & 1 # while(func.defaults()$s < 0 || func.defaults()$s > 1) { while(T) { cat("\nError: Probability out of range\nPlease re-enter a value between 0 and 1:\n") func.defaults(s = scan(what = numeric(), n = 1)) if(length(func.defaults()$s) > 0) break } } }, { # # Halfwidth precision for Heidelberger & Welch # while(T) { cat("\nEnter precision for halfwidth test (Default = 0.1):\n") func.defaults(halfwidth = scan(what = numeric(), n = 1)) if(length(func.defaults()$halfwidth) > 0) break } }, { # # Combine chains for correlation matrix # while(T) { cat("\nDo you want to combine all chains to calculate correlation matrix (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } func.defaults(combine.corr = ifelse((ans == "Y" | ans == "y"), T, F)) }, break) } } # # ----------------------------------------------------------------------- # inddat _ function(file = "bugs.ind", where = globals()$Home) { # # inddat -- function for reading BUGS index files into S-Plus # # Author: Karen Vines # Modified by Nicky Best 28/8/97 to include object Home which defines the path for the # BUGS index files (Thanks to Siem Heisterkamp for making this suggestion) # index <- scan(paste(globals()$Home, file, sep=""), list(varname = "", start = 0, finish = 0)) index } # # ----------------------------------------------------------------------- # readdat _ function(vnames, index, file = "bugs.out", start = 0, finish = 10000000, thin = 1, where = globals()$Home) { # # readdat -- function for reading BUGS output files into S-Plus # # Author: Karen Vines # Modified by Nicky Best 28/8/97 to include object Home which defines the path for the # BUGS output files (Thanks to Siem Heisterkamp for making this suggestion) # l <- 1 missed <- 1 tempa <- NULL cat("\nReading Data file...\n") temp <- scan(paste(globals()$Home, file, sep=""), list(iter = 0, value = 0)) miniter <- max(min(temp$iter), start) minadj <- (((miniter + (thin - 1)) %/% thin)) * thin maxiter <- min(max(temp$iter), finish) nx <- (maxiter %/% thin) - ((miniter + (thin - 1)) %/% thin) + 1 iter <- c(1:nx) iter <- (iter - 1) * thin + minadj nmax <- length(vnames) tempdata <- matrix(NA, nx, nmax + l) dimnames(tempdata)[2] <- list(as.character(c(1:(nmax + 1)))) tempdata[, 1] <- iter dimnames(tempdata)[[2]][1] <- c("iter") for(i in 1:nmax) { newname <- vnames[i] cat("Abstracting", newname, "... ") dimnames(tempdata)[[2]][i + l] <- newname starti <- index$start[index$varname == newname] finishi <- index$finish[index$varname == newname] if(length(starti) > 0) { tempa$value <- temp$value[starti:finishi] tempa$iter <- temp$iter[starti:finishi] useval <- tempa$value[tempa$iter %% thin == 0 & tempa$ iter <= finish & tempa$iter >= start] useiter <- tempa$iter[tempa$iter %% thin == 0 & tempa$ iter <= finish & tempa$iter >= start] tempdata[((useiter - minadj) %/% thin) + 1, i + l] <- useval cat(sum(!is.na(tempdata[, i + l])), "valid values\n") } else { cat("Not found in Index\n") } if(sum(!is.na(tempdata[, i + l])) == 0) { tempdata <- tempdata[, - (i + l)] l <- l - 1 } else { missed <- missed * is.na(tempdata[, i + l]) } } if((nmax + l) == 1) { tempdata <- NULL } else { tempdata <- tempdata[as.logical(1 - missed), ] } tempdata } # # ----------------------------------------------------------------------- # outstats _ function() { # # outstats -- Output analysis statistics for CODA system # # Author: Kate Cowles # Modified by: Nicky Best # digits <- func.defaults()$digits # # # cat("\n\nSUMMARY STATISTICS\n==================\n\n") cat("Iterations used = ", data.defaults()$start+globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="") cat("Thinning interval =", data.defaults()$thin, "\n") cat("Sample size per chain =", globals()$Niters, "\n\n") cat("Batch size for calculating Batch SE =", func.defaults()$batch.size, "\n") cat("\n\nSUMMARY STATISTICS\n==================\n\n", file = globals()$outfile, append = T) cat("Iterations used = ", data.defaults()$start+globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="", file = globals()$outfile, append = T) cat("Thinning interval =", data.defaults()$thin, "\n", file = globals()$outfile, append = T) cat("Sample size per chain =", globals()$Niters, "\n\n", file = globals()$outfile, append = T) cat("Batch size for calculating Batch SE =", func.defaults()$batch.size, "\n", file = globals()$outfile, append = T) if(globals()$Nchains>1) { if(func.defaults()$combine.stats) { cat("\nPooling over chains: ") cat("\nPooling over chains: ", file=globals()$outfile, append = T) cat(globals()$filenames[data.defaults()$chains[1]], "\n") cat(globals()$filenames[data.defaults()$chains[1]], "\n", file=globals()$outfile, append = T) for(i in 2:globals()$Nchains) { cat(" ", globals()$filenames[data.defaults()$chains[i]], "\n", sep="") cat(" ", globals()$filenames[data.defaults()$chains[i]], "\n", sep="", file=globals()$outfile, append = T) } } } cat("\n\n1. Empirical mean and standard deviation for each variable,") cat("\n plus standard error of the mean:\n\n") cat("\n\n1. Empirical mean and standard deviation for each variable,", file = globals()$outfile, append = T) cat("\n plus standard error of the mean:\n\n", file = globals()$outfile, append = T) if(func.defaults()$combine.stats) { Nchains.stats <- 1 } else Nchains.stats <- globals()$Nchains for(j in 1:Nchains.stats) { if(func.defaults()$combine.stats) { parmstats <- matrix(nrow=globals()$Nparms+3, ncol=5) for(i in 1:globals()$Nparms) { parmstats[i+3,1] <- signif(mean(globals()$parmsmat[,i,]), digits) parmstats[i+3,2] <- signif(sqrt(var(as.vector( globals()$parmsmat[,i,]))), digits) parmstats[i+3,3] <- sqrt(var(as.vector( globals()$parmsmat[,i,]))/length(globals()$parmsmat[,i,])) parmstats[i+3,4:5] <- calcbatch2(globals()$parmsmat[,i,], func.defaults()$batch.size)[2:3] } parmstats[-1:-3,] <- format(signif(as.numeric(parmstats[-1:-3,]), digits=digits)) parmstats[1:3, 1:5] <- c("Mean", "====", "", "SD", "==", "", "Naive SE", "========", "", "Batch SE", "========", "", "Lag-1 batch autocorr", "====================", "") } else { parmstats <- matrix(nrow=globals()$Nparms+3, ncol=6) for(i in 1:globals()$Nparms) { parmstats[i+3,1] <- signif(mean(globals()$parmsmat[,i,j]), digits) parmstats[i+3,2] <- signif(sqrt(var(as.vector( globals()$parmsmat[,i,j]))), digits) parmstats[i+3,3] <- sqrt(var(as.vector( globals()$parmsmat[,i,j]))/length(globals()$parmsmat[,i,j])) parmstats[i+3,4] <- geweke.nse(globals()$parmsmat[,i,j]) parmstats[i+3,5:6] <- calcbatch2(globals()$parmsmat[,i,j], func.defaults()$batch.size)[2:3] } parmstats[-1:-3,] <- format(signif(as.numeric(parmstats[-1:-3,]), digits=digits)) parmstats[1:3, 1:6] <- c("Mean", "====", "", "SD", "==", "", "Naive SE", "========", "","Time-series SE", "==============", "", "Batch SE", "========", "", "Lag-1 batch autocorr", "====================", "") cat("\nChain:", globals()$filenames[data.defaults()$chains[j]], "\n") cat(paste(rep("=", 7 + nchar( globals()$filenames[data.defaults()$chains[j]])), collapse=""), "\n\n") cat("\nChain:", globals()$filenames[data.defaults()$chains[j]], "\n", file=globals()$outfile, append = T) cat(paste(rep("=", 7 + nchar( globals()$filenames[data.defaults()$chains[j]])), collapse=""), "\n\n", file=globals()$outfile, append = T) } # # Calculate how many columns will fit on page - if more columns # are needed than will fit, then spread over 2 or more tables # cwidth <- 80 - (max(nchar(globals()$Labels)) + 6) colsize <- apply(nchar(parmstats), 2, max) + 5 tmp <- cumsum(colsize) nr <- list() nr[[1]] <- c(1:length(colsize))[tmp <= cwidth] i<-1 while(nr[[i]][length(nr[[i]])] != length(colsize)) { tmp <- tmp + (cwidth*i - tmp[nr[[i]][length(nr[[i]])]]) nr[[i+1]] <- c(1:length(colsize))[tmp > (cwidth*i) & tmp <= (cwidth*(i+1))] i <- i+1 } ntabs <- length(nr) for(n in 1:ntabs) { stats.table<- matrix(globals()$parmnames, ncol=1) stats.table <- cbind(stats.table, matrix( pretty.print(parmstats[,nr[[n]]], print.it=F, hsep=" ", vsep="",csep=" "), ncol=1)) stats.table <- pretty.print(stats.table, hsep=" | ", vsep="-", csep="-+-", print.it=F) cat(stats.table) cat(stats.table, file=globals()$outfile, append = T) } } # # # cat("\n\n2. Quantiles for each variable:\n\n") cat("\n\n2. Quantiles for each variable:\n\n", file = globals()$outfile, append = T) quant <- func.defaults()$quantiles parmstats <- matrix(nrow=globals()$Nparms+3, ncol=length(quant)) for(k in 1:Nchains.stats) { for(i in 1:globals()$Nparms) { if(func.defaults()$combine.stats) { parmstats[i+3,] <- quantile(as.vector(globals()$parmsmat[,i,]), quant) } else { parmstats[i+3,] <- quantile(as.vector(globals()$parmsmat[,i,k]), quant) } } quant.labels <- paste(100*quant, "%", sep="") for(j in 1:length(quant)) { parmstats[1:3, j] <- c(quant.labels[j], paste(rep("=", nchar( quant.labels[j])), collapse=""), "") } parmstats[-1:-3,] <- format(signif(as.numeric(parmstats[-1:-3,]), digits=digits)) if(!func.defaults()$combine.stats) { cat("\nChain:", globals()$filenames[data.defaults()$chains[k]], "\n") cat(paste(rep("=", 7+nchar(globals()$filenames[data.defaults()$chains[k]])), collapse=""), "\n\n") cat("\nChain:", globals()$filenames[data.defaults()$chains[k]], "\n", file=globals()$outfile, append=T) cat(paste(rep("=", 7+nchar(globals()$filenames[data.defaults()$chains[k]])), collapse=""), "\n\n", file=globals()$outfile, append=T) } # # Calculate how many columns will fit on page - if more columns # are needed than will fit, then spread over 2 or more tables # cwidth <- 80 - (max(nchar(globals()$Labels)) + 6) colsize <- apply(nchar(parmstats), 2, max) + 5 tmp <- cumsum(colsize) nr <- list() nr[[1]] <- c(1:length(colsize))[tmp <= cwidth] i<-1 while(nr[[i]][length(nr[[i]])] != length(colsize)) { tmp <- tmp + (cwidth*i - tmp[nr[[i]][length(nr[[i]])]]) nr[[i+1]] <- c(1:length(colsize))[tmp > (cwidth*i) & tmp <= (cwidth*(i+1))] i <- i+1 } ntabs <- length(nr) for(n in 1:ntabs) { stats.table <- matrix(globals()$parmnames, ncol=1) stats.table <- cbind(stats.table, matrix(pretty.print(parmstats[,nr[[n]]], print.it=F, hsep=" ", vsep="",csep=" "), ncol=1)) stats.table <-pretty.print(stats.table, hsep=" | ", vsep="-", csep="-+-", print.it=F) cat(stats.table) cat(stats.table, file=globals()$outfile, append = T) } } } # # ----------------------------------------------------------------------- # calcbatch2 _ function(a, batchsiz) { # # calcbatch2 -- standard errors by batch means method # # Author: Kate Cowles # l <- length(a) k <- batchsiz m <- floor(l/k) # number of batches meanvect <- numeric() for(i in 0:(m - 1)) { meanvect[i + 1] <- mean(a[(k * i + 1):(k * (i + 1))]) } myvar <- numeric() myvar[1] <- mean(meanvect) myvar[2] <- sqrt(var(meanvect)/m) myvar[3] <- acf(meanvect, 1, "corr", F)$acf[2] myvar[4] <- k names(myvar) <- c("Mean", "SEM", "Lag 1 autocorr between batches", "Batch size") myvar } # # ----------------------------------------------------------------------- # outplots _ function() { # # outplots -- Plots for CODA system # # Author: Kate Cowles # Modified by: Nicky Best # open.graphics() call.plots(mtitle = globals()$mtitle, densplot = func.defaults()$densplot, trace = func.defaults()$trace, use.lowess = func.defaults()$lowess, which.pos = func.defaults()$pos.parmsmat, which.prop = func.defaults()$prop.parmsmat, sepplot = func.defaults()$sepplot) } # # ----------------------------------------------------------------------- # setup.plots _ function(Nparms.tmp, nplots = 1, sepplot = F, which.plot, args.list, one.page = func.defaults()$onepage, extra.title = NULL) { # # setup.plots -- sets up graphics window and options for plotting raw BUGS # output, Gelman & Rubin diagnostics or autocorrelations # # Author: Nicky Best # plots <- T multiplot <- F if(one.page) { Nparms.setup <- 1 } else Nparms.setup <- Nparms.tmp if(!func.defaults()$user.layout) { # # 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(sepplot && globals()$Nchains > 1) { if(nplots == 1) { mrows <- ifelse((Nparms.setup == 1 && globals()$Nchains == 2), 1, ifelse((Nparms.setup == 2 && globals()$Nchains == 2) || (Nparms.setup == 1 && any(globals()$Nchains == c(3, 4))) || (any(Nparms.setup == c(2,4)) && globals()$Nchains == 3), 2, ifelse((Nparms.setup == 4 && globals()$Nchains == 2) || (globals()$Nchains == 4 && Nparms.setup > 1), 4, 3))) mcols <- ifelse((globals()$Nchains == 3 && any(Nparms.setup == c(2, 4))) || ((Nparms.setup == 3 || Nparms.setup >= 5) && globals()$Nchains == 3) || any(globals()$Nchains == c(7:9)) || globals()$Nchains >= 13, 3, 2) } else { mrows <- ifelse(globals()$Nchains <= 4, globals()$Nchains, ifelse(any(globals()$Nchains == c(5,6,10:12)), 3, 4)) mcols <- 2 } } else { if(nplots == 1) { mrows <- ifelse(Nparms.setup <= 2, 1, ifelse(Nparms.setup <= 4, 2, 3)) mcols <- ifelse(Nparms.setup <= 9, Nparms.setup %/% mrows + ifelse(Nparms.setup %% mrows == 0, 0, 1), ifelse(any(Nparms.setup == c(10:12)), 2, 3)) } else { mrows <- ifelse(Nparms.setup <= 4, Nparms.setup, ifelse(any(Nparms.setup == c(5,6,10:12)), 3, 4)) mcols <- 2 } } } else { mrows <- func.defaults()$mrows mcols <- func.defaults()$mcols } if(which.plot[1]=="parmcorr.plot") { args.list$mrows <- mrows } 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 while(plots) { if(globals()$which.dev < 4) { # # UNIX options: # ============= # if(globals()$which.dev == 1 | globals()$which.dev ==2) { # # Openlook & Motif # par(mfrow = c(mrows, mcols), oma = c(1.1, 0, 3.8, 0), mgp = c(3.9, .9, 0)) par(mar = c(6.1, 6.5, 5.1, 3.1)) } else { # # X11 # par(mfrow = c(mrows, mcols), oma = c(1.3, 0, 4.2, 0), mgp = c(4.2, .9, 0)) par(mar = c(6.7, 7, 5.7, 3.5)) } args.list$scale <- switch(plot.dim, 1.1, 1.0, 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 ) } else { # # PC options: # =========== # par(mfrow = c(mrows, mcols), oma = c(1, 0, 3.8, 0), mgp = c(1.5, 0.4, 0)) par(mar = c(3.25, 3.5, 1.5, 2)) args.list$scale <- switch(plot.dim, 0.9, 0.8, 0.7, 0.65, 0.6, 0.5, 0.425, 0.375, 0.2, 0.2) args.list$title.scale <- switch(plot.dim, 0.65, 0.6, 0.55, 0.5, 0.45, 0.4, 0.25, 0.15, 0.15, 0.1) } if(plot.dim<3 & func.defaults()$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)) } } if(nplots==2) { p1 <- 1 p2 <- 2 } else { if(func.defaults()$trace) { p1 <- 1 p2 <- 1 } else { p1 <- 2 p2 <- 2 } } for(j in 1:Nparms.tmp) { # # Select jth variable for plotting # args.list$this.parm <- j for(k in p1:p2) { if(sepplot && globals()$Nchains > 1) { for(i in 1:globals()$Nchains) { args.list$this.chain <- i do.call(eval(which.plot[k]), args.list) # # If this is last plot on page, add title: # if((par()$mfg[1] == par()$mfg[3] && par()$mfg[2] == par()$mfg[4]) || (k==p2 && i==globals()$Nchains && j==Nparms.tmp) || (k==p2 && i==globals()$Nchains && one.page)) { if(length(globals()$mtitle) > 0) { if(func.defaults()$ps.plot < 3) { mtext(side = 3, line = 1.35, cex = 1.1, outer = TRUE, globals()$mtitle) } else { mtext(side = 3, line = 1.7, cex = 0.8, outer = TRUE, globals()$mtitle) } if(length(extra.title) > 0) { if(func.defaults()$ps.plot < 3) { mtext(side = 3, line = -0.35, cex=0.85, outer = TRUE, extra.title) } else { mtext(side = 3, line = -0.4, cex=0.65, outer = TRUE, extra.title) } } } else { if(length(extra.title) > 0) { if(func.defaults()$ps.plot < 3) { mtext(side = 3, line = 1.35, cex = 1.1, outer = TRUE, extra.title) } else { mtext(side = 3, line = 1.35, cex = 0.9, outer = TRUE, extra.title) } } } # # If plots span more than 1 page, insert a # pause before displaying next page of plots: # if(j < Nparms.tmp || (j == Nparms.tmp && i < globals()$Nchains) || (Nparms.setup == 1 && i < globals()$Nchains)) { multiplot <- T pause() } # # Ensure that next plot is started on new page: # if(globals()$which.dev < 4) { # # UNIX options: # ============= # if(globals()$which.dev == 1 | globals()$which.dev ==2) { # # Openlook & Motif # par(mfrow = c(mrows, mcols), oma = c(1.1, 0, 3.8, 0), mgp = c(3.9, .9, 0)) par(mar = c(6.1, 6.5, 5.1, 3.1)) } else { # # X11 # par(mfrow = c(mrows, mcols), oma = c(1.3, 0, 4.2, 0), mgp = c(4.1, .9, 0)) par(mar = c(6.7, 7, 5.7, 3.5)) } } else { # # PC options: # =========== # par(mfrow = c(mrows, mcols), oma = c(1, 0, 3.8, 0), mgp = c(1.5, 0.4, 0)) par(mar = c(3.25, 3.5, 1.5, 2)) } } } } else { args.list$this.chain <- 1 do.call(eval(which.plot[k]), args.list) # # If this is last plot on page, add title: # if((par()$mfg[1] == par()$mfg[3] && par()$mfg[2] == par()$mfg[4]) || (k==p2 && j==Nparms.tmp) || (k==p2 && one.page)) { if(length(globals()$mtitle) > 0) { if(func.defaults()$ps.plot < 3) { mtext(side = 3, line = 1.35, cex = 1.1, outer = TRUE, globals()$mtitle) } else { mtext(side = 3, line = 1.7, cex = 0.8, outer = TRUE, globals()$mtitle) } if(length(extra.title) > 0) { if(func.defaults()$ps.plot < 3) { mtext(side = 3, line = -0.35, cex=0.85, outer = TRUE, extra.title) } else { mtext(side = 3, line = -0.4, cex=0.65, outer = TRUE, extra.title) } } } else { if(length(extra.title) > 0) { if(func.defaults()$ps.plot < 3) { mtext(side = 3, line = 1.35, cex = 1.1, outer = TRUE, extra.title) } else { mtext(side = 3, line = 1.35, cex = 0.9, outer = TRUE, extra.title) } } } # # If plots span more than 1 page, insert a # pause before displaying next page of plots # if(j < Nparms.tmp) { multiplot <- T pause() } } } } } # # After last plot is completed, query whether plots should # be saved as a postscript file # plot.to.file() # # If multiple pages of plots have been displayed, query whether # user wants to return to first plot again (e.g. to print it), or # continue with the Output Analysis Menu # if(multiplot) { while(T) { cat("\nFinished displaying requested plots\nDo you want to display the first page of plots again (y/n)?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } plots <- ifelse((ans == "Y" | ans == "y"), T, F) } else { plots <- F } } } # # ----------------------------------------------------------------------- # call.plots _ function(mtitle, densplot = T, trace = T, use.lowess = T, which.pos = 0, which.prop = 0, sepplot = F) { # # call.plots -- function to set up required arguments for calling draw.trace() # and draw.dens() # # Author: Nicky Best # v <- globals()$Labels[data.defaults()$parms] nplots <- NULL as.integer(nplots) nplots <- densplot + trace if(nplots == 0) stop("No plots requested") if(length(v) == 0) stop("No variables specified") # # Identify positive-valued parameters # pos <- rep(FALSE, length(v)) if(length(which.pos) > 1 || (length(which.pos) == 1 && which.pos > 0)) { for(k in 1:length(which.pos)) { if(which.pos[k] > length(v)) stop("Positive variable indicator out of range" ) pos[which.pos[k]] <- TRUE } } # # Identify parameters taking values on (0, 1) # prop <- rep(FALSE, length(v)) if(length(which.prop) > 1 || (length(which.prop) == 1 && which.prop > 0)) { for(k in 1:length(which.prop)) { if(which.prop[k] > length(v)) stop("(0, 1) variable indicator out of range" ) prop[which.prop[k]] <- TRUE } } args.list <- list() args.list[[1]] <- sepplot args.list[[2]] <- numeric() args.list[[3]] <- numeric() args.list[[4]] <- func.defaults()$densplot args.list[[5]] <- func.defaults()$trace args.list[[6]] <- func.defaults()$lowess args.list[[7]] <- pos args.list[[8]] <- prop args.list[[9]] <- numeric() args.list[[10]] <- numeric() args.list[[11]] <- numeric() args.list[[12]] <- numeric() args.list[[13]] <- numeric() args.list[[14]] <- numeric() args.list[[15]] <- numeric() args.list[[16]] <- numeric() args.list[[17]] <- numeric() args.list[[18]] <- numeric() if(trace) { # # If multiple chains, set scale to be same for all plots # of each variable # args.list[[11]] <- apply(globals()$parmsmat, 2, min) args.list[[12]] <- apply(globals()$parmsmat, 2, max) if(data.defaults()$all.trace) { # need to calculate min and max of ALL iterations for trace plots args.list[[13]] <- apply(globals()$out[[data.defaults()$chains[1]]], 2, min) args.list[[14]] <- apply(globals()$out[[data.defaults()$chains[1]]], 2, max) # need to calculate max of 2nd half iterations for legend # on trace plots n2 <- dim(globals()$out[[1]])[1] n1 <- floor(n2/2) args.list[[15]] <- apply(globals()$out[[data.defaults()$chains[1]]][n1:n2,], 2, max) for(j in 1:globals()$Nchains) { tmp1 <- apply(globals()$out[[data.defaults()$chains[j]]], 2, min) tmp2 <- apply(globals()$out[[data.defaults()$chains[j]]], 2, max) for(i in 1:length(tmp1)) { args.list[[13]][i] <- min(args.list[[13]][i], tmp1[i]) args.list[[14]][i] <- max(args.list[[14]][i], tmp2[i]) } tmp1 <- apply(globals()$out[[data.defaults()$chains[j]]][n1:n2,], 2, max) for(i in 1:length(tmp1)) { args.list[[15]][i] <- max(args.list[[15]][i], tmp1[i]) } } args.list[[13]] <- args.list[[13]][-1] args.list[[14]] <- args.list[[14]][-1] args.list[[13]] <- args.list[[13]][data.defaults()$parms] args.list[[14]] <- args.list[[14]][data.defaults()$parms] args.list[[15]] <- args.list[[15]][-1] args.list[[15]] <- args.list[[15]][data.defaults()$parms] } else { args.list[[13]] <- args.list[[11]] args.list[[14]] <- args.list[[12]] if(globals()$Nparms==1 & globals()$Nchains==1) { args.list[[15]] <- max(globals()$parmsmat[floor(globals()$Niters/2):globals()$Niters,,]) } else{ args.list[[15]] <- apply(globals()$parmsmat[floor(globals()$Niters/2):globals()$Niters,,], 2, max) } } } if(densplot) { # # Now calculate scale for density plots # if(sepplot) { for(j in 1:globals()$Nchains) { tmp1 <- apply(as.matrix(globals()$parmsmat[, , j]), 2, min) tmp2 <- apply(as.matrix(globals()$parmsmat[, , j]), 2, max) tmp3 <- apply(as.matrix(globals()$parmsmat[, , j]), 2, func.defaults()$bandwidth) tmp3[tmp3==0] <- 0.1 tmp4<-numeric() for(i in 1:globals()$Nparms) { if(pos[i]) { tmpy <- c(globals()$parmsmat[, i, j], -globals()$parmsmat[, i, j]) } else { if(prop[i]) { tmpy <- c(globals()$parmsmat[, i, j], -globals()$parmsmat[, i, j], 2-globals()$parmsmat[, i, j]) } else tmpy <- globals()$parmsmat[, i, j] } tmp4[i] <- max((density(tmpy, width=tmp3[i], n=200))$y) if(pos[i]) tmp4[i] <- 2*tmp4[i] if(prop[i]) tmp4[i] <- 3*tmp4[i] } if(j==1) { args.list[[16]] <- tmp1 - 0.75*tmp3 args.list[[17]] <- tmp2 + 0.75*tmp3 args.list[[18]] <- tmp4 } else { tmp11 <- rbind(args.list[[16]], (tmp1 - 0.75*tmp3)) tmp22 <- rbind(args.list[[17]], (tmp2 + 0.75*tmp3)) args.list[[16]] <- apply(tmp11, 2, min) args.list[[17]] <- apply(tmp22, 2, max) tmp44 <- rbind(args.list[[18]], tmp4) args.list[[18]] <- apply(tmp44, 2, max) } } } else { tmp1 <- apply(globals()$parmsmat, 2, min) tmp2 <- apply(globals()$parmsmat, 2, max) tmp3 <- apply(globals()$parmsmat, 2, func.defaults()$bandwidth) tmp3[tmp3==0] <- 0.1 args.list[[16]] <- tmp1 - 0.75*tmp3 args.list[[17]] <- tmp2 + 0.75*tmp3 args.list[[18]] <- numeric() for(i in 1:globals()$Nparms) { if(pos[i]) { tmpy <- c(globals()$parmsmat[, i, ], -globals()$parmsmat[, i, ]) } else { if(prop[i]) { tmpy <- c(globals()$parmsmat[, i, ], -globals()$parmsmat[, i, ], 2-globals()$parmsmat[, i, ]) } else tmpy <- globals()$parmsmat[, i, ] } args.list[[18]][i] <- max((density(tmpy, width=tmp3[i], n=200))$y) if(pos[i]) args.list[[18]][i] <- 2*args.list[[18]][i] if(prop[i]) args.list[[18]][i] <- 3*args.list[[18]][i] } } } names(args.list) <- c("sepplot", "this.parm", "this.chain", "densplot", "trace", "use.lowess", "pos", "prop", "scale", "title.scale", "minlim", "maxlim", "minlim2", "maxlim2","maxlim3","densxmin", "densxmax", "densymax") setup.plots(globals()$Nparms, nplots, sepplot, which.plot = c("draw.trace", "draw.dens"), args.list) } # # ----------------------------------------------------------------------- # draw.trace _ function(sepplot=F, this.parm=1, this.chain=1, densplot=T, trace=T, use.lowess=T, pos, prop, scale=1, title.scale=1, minlim, maxlim, minlim2, maxlim2, maxlim3, densxmin, densxmax, densymax) { # # draw.trace -- function for plotting traces of BUGS output # # Author: Karen Vines # Modified by: Nicky Best # if(trace) { v <- globals()$Labels[data.defaults()$parms] cn <- globals()$filenames[data.defaults()$chains] leg.scale <- scale*0.6 if(func.defaults()$ps.plot == 3) { leg.scale <- scale*0.6 title.scale <- title.scale*.7 scale <- scale*.6 } if(data.defaults()$all.trace) { # need to use all iterations for trace plot, so take from out # rather than parmsmat, as latter may not contain all samples if(sepplot) { ytemp <- globals()$out[[data.defaults()$chains[this.chain]]][, data.defaults()$parms[this.parm]+1] xtemp <- globals()$out[[data.defaults()$chains[this.chain]]][,1] } else { ytemp <- numeric() for(i in 1:globals()$Nchains) { ytemp <- c(ytemp, globals()$out[[data.defaults()$chains[i]]][, data.defaults()$parms[this.parm]+1]) } xtemp <- rep(globals()$out[[data.defaults()$chains[this.chain]]][,1], length(data.defaults()$chains)) } } else { if(sepplot) { ytemp <- globals()$parmsmat[, this.parm, this.chain] xtemp <- globals()$parmsmat.iter } else { ytemp <- globals()$parmsmat[,this.parm,] xtemp <- rep(globals()$parmsmat.iter, length(data.defaults()$chains)) } } miss <- !is.na(ytemp) y <- ytemp[miss] x <- xtemp[miss] n <- length(y) if(n > 0) { # determine y limits for plot depending on max(y) # over all chains, and number of lines of legend maxlim3[this.parm] <- maxlim3[this.parm] + (0.45 + 0.15*globals()$Nchains)* (maxlim3[this.parm] - minlim2[this.parm]) maxlim3[this.parm] <- max(maxlim2[this.parm], maxlim3[this.parm]) if(sepplot) { tmp <- c(minlim2[this.parm],maxlim2[this.parm]) lim <- range(pretty(tmp, 2)) plot(x, y, type = "l", xlab = "iteration", ylab = v[this.parm], cex = scale, ylim=lim, axes=F) axis(1, at=pretty(x,3), labels=my.format(pretty(x,3), scientific=c(-func.defaults()$scientific, func.defaults()$scientific)), cex = scale) par(crt=90, srt=90) axis(2, at=pretty(tmp, 2), labels=my.format(pretty(tmp, 2), scientific=c(-func.defaults()$scientific, func.defaults()$scientific)), cex = scale) par(crt=0, srt=0) box() mtext(paste("Trace of ", cn[this.chain], sep = ""), cex = title.scale*1.3, side=3, line=2.5) mtext(paste("(",format(n), " values)", sep=""), cex = title.scale*1.3, side=3, line=.65) if(use.lowess) { lines(lowess(x, y)) } } else { ymax <- ifelse(globals()$Nchains > 1, maxlim3[this.parm], max(y)) if(data.defaults()$all.trace) { nn <- dim(globals()$out[[1]])[1] } else nn <- dim(globals()$parmsmat)[1] tmp <- c(minlim2[this.parm], ymax) lim <- range(pretty(tmp, 2)) plot(x, y, type = "n", xlab = "iteration", ylab = v[this.parm], cex = scale, ylim=lim, axes=F) axis(1, at=pretty(x,3), labels=my.format(pretty(x,3), scientific=c(-func.defaults()$scientific, func.defaults()$scientific)), cex = scale) par(crt=90, srt=90) axis(2, at=pretty(tmp, 2), labels=my.format(pretty(tmp, 2), scientific=c(-func.defaults()$scientific, func.defaults()$scientific)),cex = scale) par(crt=0, srt=0) box() for(i in 1:length(data.defaults()$chains)) { lines(x[(nn*(i-1)+1):(nn*i)], y[(nn*(i-1)+1):(nn*i)], lty=i) if(use.lowess) { lines(lowess( x[(nn*(i-1)+1):(nn*i)], y[(nn*(i-1)+1):(nn*i)]), lty=i) } } if(globals()$Nchains > 1) { cnlab <- character(globals()$Nchains) ltlab <- numeric(globals()$Nchains) for(i in 1:globals()$Nchains) { cnlab[i] <- cn[i] ltlab[i] <- i } loc <- c(0.7, 0.48, 0.26) nc <- globals()$Nchains %/% 2 + globals()$Nchains %% 2 nc <- ifelse(nc > 3, 3, nc) legend(min(x) + (max(x) - min(x)) * loc[nc], lim[2] + 0.02 * (lim[2] - lim[1]), cnlab, lty = ltlab, cex = 0.75*leg.scale, bty = "n", ncol = nc) } mtext(paste("Trace of ", v[this.parm], sep = ""), cex = title.scale*1.3, side=3, line=2.5) mtext(paste("(",format(n/ length(data.defaults()$chains)), " values per trace)", sep=""), cex = title.scale*1.3, side=3, line=.65) } } else if(n == 0) { cat(variable, "has no valid values \n") } } } # # ----------------------------------------------------------------------- # draw.dens _ function(sepplot=F, this.parm=1, this.chain=1, densplot=T, trace=T, use.lowess=T, pos, prop, scale=1, title.scale=1, minlim, maxlim, minlim2, maxlim2, maxlim3, densxmin, densxmax, densymax) { # # draw.dens -- function for plotting density estiamtes # # Author: Karen Vines # Modified by: Nicky Best # if(densplot) { v <- globals()$Labels[data.defaults()$parms] cn <- globals()$filenames[data.defaults()$chains] leg.scale <- scale*0.6 if(func.defaults()$ps.plot == 3) { leg.scale <- scale*0.6 title.scale <- title.scale*.7 scale <- scale*.6 } if(sepplot) { ytemp <- globals()$parmsmat[, this.parm, this.chain] xtemp <- globals()$parmsmat.iter } else { ytemp <- globals()$parmsmat[,this.parm,] xtemp <- rep(globals()$parmsmat.iter, length(data.defaults()$chains)) } miss <- !is.na(ytemp) y <- ytemp[miss] yy <- y x <- xtemp[miss] n <- length(y) if(n > 0) { xlimits <- c(densxmin[this.parm], densxmax[this.parm]) ylimits <- c(0, densymax[this.parm]) if(pos[this.parm]) { xlimits <- c(0, densxmax[this.parm]) } if(prop[this.parm]) { xlimits <- c(0, 1) } ww <- func.defaults()$bandwidth(y) # # Occassionally, variables corresponding to # proportions/probabilities are repeatedly # sampled as 0 or 1 (i.e. every value is the same). # => ww = 0, so need to replace this with a small # positive value, say 0.1 # if(ww == 0) { ww <- 0.1 } if(pos[this.parm]) { # Reflect all points in boundary y=0 y <- c(y, -y) } if(prop[this.parm]) { # Reflect all points in boundaries y=0 & y=1 y <- c(y, -y, 2-y) } dens <- density(y, width = ww, n=200) if(pos[this.parm]) { # f(y) = 2*dens$y for y>=0; 0 otherwise dens$y[dens$x>=0] <- 2*dens$y[dens$x>=0] dens$y[dens$x<0] <- 0 dens$x[dens$x<0] <- 0 } if(prop[this.parm]) { # f(y) = 3*dens$y for y>=0 & y<=1; 0 otherwise 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 } plot(dens, type = "l", xlab = v[this.parm], ylab = "", cex = scale, xlim=xlimits, ylim=ylimits, axes=F) axis(1, at=pretty(xlimits,3), labels=my.format(pretty(xlimits,3), scientific=c(-func.defaults()$scientific, func.defaults()$scientific)), cex = scale) par(crt=90, srt=90) axis(2, at=pretty(ylimits, 2), labels=my.format(pretty(ylimits, 2), scientific=c(-func.defaults()$scientific, func.defaults()$scientific)),cex = scale) par(crt=0, srt=0) box() if(sepplot) { mtext(paste("Kernel density from ", cn[this.chain], sep = ""), cex = title.scale*1.3, side=3, line=2.5) mtext(paste("(",format(n), " values)", sep=""), cex = title.scale*1.3, side=3, line=.65) } else { mtext(paste("Kernel density for ", v[this.parm], sep = ""), cex = title.scale*1.3, side=3, line=2.5) mtext(paste("(",format(n), " values)", sep=""), cex = title.scale*1.3, side=3, line=.65) } points(yy, yy * 0) } else if(n == 0) { cat(variable, "has no valid values \n") } } } # # ----------------------------------------------------------------------- # randl2 _ function() { # # randl2 -- CODA system module to write files for calculating Raftery & Lewis's # diagnostic # # Author: Kate Cowles # Modified by: Nicky Best # { cat("\nRAFTERY AND LEWIS CONVERGENCE DIAGNOSTIC:") cat("\n=========================================\n\n") cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="") cat("Thinning interval =", data.defaults()$thin, "\n") cat("Sample size per chain =", globals()$Niters, "\n") cat("\nQuantile =", func.defaults()$q) cat("\nAccuracy = +/-", func.defaults()$r) cat("\nProbability =", func.defaults()$s, "\n") cat("\nRAFTERY AND LEWIS CONVERGENCE DIAGNOSTIC:\n", file = globals()$outfile, append = T) cat("\n=========================================\n\n", file = globals()$outfile, append = T) cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="", file = globals()$outfile, append = T) cat("Thinning interval =", data.defaults()$thin, "\n", file = globals()$outfile, append = T) cat("Sample size per chain =", globals()$Niters, "\n", file = globals()$outfile, append = T) cat("\nQuantile =", func.defaults()$q, file = globals()$outfile, append = T) cat("\nAccuracy = +/-", func.defaults()$r, file = globals()$outfile, append = T) cat("\nProbability =", func.defaults()$s, "\n", file = globals()$outfile, append = T) for(j in 1:globals()$Nchains) { RL.mat <- gibbsit2(matrix(globals()$parmsmat[, , j], nrow=dim(globals()$parmsmat)[1], ncol=dim(globals()$parmsmat)[2]), globals()$Labels[data.defaults()$parms], q=func.defaults()$q, r=func.defaults()$r, s=func.defaults()$s) cat("\n\nChain: ", globals()$filenames[data.defaults()$chains[j]], "\n", paste(c("=======", rep("=", nchar(globals()$filenames[data.defaults()$chains[j]]))), collapse=""), "\n\n", sep="") cat("\n\nChain: ", globals()$filenames[data.defaults()$chains[j]], "\n", paste(c("=======", rep("=", nchar(globals()$filenames[data.defaults()$chains[j]]))), collapse=""), "\n\n", sep="", file=globals()$outfile, append=T) if(RL.mat[1,1] == "Error") { cat("\n******* Error: *******\nNmin = ", RL.mat[2,1], "; Available chain length is ", globals()$Niters, "\nRe-run BUGS program for at least ", RL.mat[2,1], " iterations\nOR reduce accuracy, probability and/or quantile to be estimated\n\n", sep="") cat("\n******* Error: *******\nNmin = ", RL.mat[2,1], "; Available chain length is ", globals()$Niters, "\nRe-run BUGS program for at least ", RL.mat[2,1], " iterations\nOR reduce accuracy, probability and/or quantile to be estimated\n\n", sep="", file=globals()$outfile, append=T) } else { nmin <- RL.mat[1,4] RL.mat <- rbind(matrix(c("Thin", "Burn-in ", "Total", "Lower bound ", "Dependence", "(k)", " (M)", "(N)", " (Nmin)", "factor (I)","====", "=======", "=====", "===========", "==========", "", "", "", "", ""), byrow=T, ncol=5), RL.mat) RL.table <- pretty.print(matrix(c(globals()$parmnames2, pretty.print(RL.mat, print.it=F, hsep=" ", vsep="", csep=" ")), ncol=2), hsep=" | ", vsep="-", csep="-+-", print.it=F) cat(RL.table) cat(RL.table, file=globals()$outfile, append=T) } } } } # # ----------------------------------------------------------------------- # gibbsit2 _ function(data, vnames, q = 0.025, r = 0.005, s = 0.95, delta = 0.001) { # # gibbsit2 -- Code for Raftery And Lewis Diagnostic # Estimates the theoretical minimum number of iterations (nmin) # plus the actual number of initial iterations to be discarded # (nburn) and the number of additional iterations to keep (nkeep), # based on the output from BUGS run, in order to estimate the # cumulative distribution function of the qth quantile of each # variable to within +/- r with probability s. The quantity # # I = (nburn + nkeep) / nmin # # measures the increase in number of iterations needed due to # dependence in the sampled sequence. Values much greater than # 1.0 (~ >5.0) indicate high within-chain correlations and # probable convergence problems # # Author: Karen Vines # Modified by: Nicky Best # resmatrix <- matrix(nrow=globals()$Nparms, ncol=5) miniter <- globals()$parmsmat.iter[1] maxiter <- globals()$parmsmat.iter[nrow(data)] # # nmin = minimun number of iterations required to estimate q to within # +/- r with probability s if successive values of Z_i were independent # where Z_i = d(U_i <= u); U_i is the value the variable of interest at # iteration i; u is the sampled value of the quantile q for this variable; # and d is an indicator function) # phi <- qnorm(0.5 * (1 + s)) nmin <- as.integer(ceiling((q * (1 - q) * phi^2)/r^2)) iteracnt <- nrow(data) if(nmin > globals()$Niters) { # # Error: nmin > number of iterations # => don't compute remaining diagnostics # resmatrix <- matrix(c("Error", nmin), ncol=1) } else { # # Initializing... # for(i in 1:globals()$Nparms) { # # u = quant # Z_i = 1 if data[, i] <= quant; 0 otherwise # quant <- quantile(data[, i], probs = q) dichot <- as.integer(data[,i] <= quant) kwork <- 0 bic <- 1 # # First need to find the thinning parameter kthin # while(bic >= 0) { kwork <- as.integer(kwork + 1) testres <- dichot[seq(1, iteracnt, by = kwork)] newdim <- length(testres) tttemp <- as.integer(testres[1:(newdim - 2)] + testres[ 2:(newdim - 1)] * 2 + testres[3:newdim] * 4) testtran <- array(as.integer(0), dim=c(2,2,2)) testtran[1,1,1] <- sum(as.integer(tttemp == 0)) testtran[2,1,1] <- sum(as.integer(tttemp == 1)) testtran[1,2,1] <- sum(as.integer(tttemp == 2)) testtran[2,2,1] <- sum(as.integer(tttemp == 3)) testtran[1,1,2] <- sum(as.integer(tttemp == 4)) testtran[2,1,2] <- sum(as.integer(tttemp == 5)) testtran[1,2,2] <- sum(as.integer(tttemp == 6)) testtran[2,2,2] <- sum(as.integer(tttemp == 7)) 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 * data.defaults()$thin) # # then need to find length of burn-in and for required precision # fttemp <- as.integer(testres[1:(newdim - 1)] + testres[2: newdim] * 2) fttable <- table(fttemp) finaltran <- matrix(as.integer(0), nrow = 2, ncol = 2) if(!is.na(fttable["0"])) finaltran[1,1] <- fttable["0"] if(!is.na(fttable["1"])) finaltran[2,1] <- fttable["1"] if(!is.na(fttable["2"])) finaltran[1,2] <- fttable["2"] if(!is.na(fttable["3"])) finaltran[2,2] <- fttable["3"] alpha <- finaltran[1,2]/(finaltran[1,1] + finaltran[1,2]) beta <- finaltran[2,1]/(finaltran[2,1] + finaltran[2,2]) tempburn <- log((delta * (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] <- kthin resmatrix[i, 2] <- nburn resmatrix[i, 3] <- nkeep+nburn resmatrix[i, 4] <- nmin resmatrix[i, 5] <- signif(iratio, digits=3) } } resmatrix } # # ----------------------------------------------------------------------- # geweke _ function() { # # geweke -- Calculate Geweke convergence diagnostic for CODA system # # Author: Kate Cowles # Modified by: Nicky Best # cat("\nGEWEKE CONVERGENCE DIAGNOSTIC (Z-score):") cat("\n========================================\n\n") cat("Iterations used = ", data.defaults()$start+globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="") cat("Thinning interval =", data.defaults()$thin, "\n") cat("Sample size per chain =", globals()$Niters, "\n") cat("\nFraction in 1st window =", func.defaults()$frac1) cat("\nFraction in 2nd window =", func.defaults()$frac2, "\n\n") cat("\nGEWEKE CONVERGENCE DIAGNOSTIC (Z-score):", file = globals()$outfile, append = T) cat("\n========================================\n\n", file = globals()$outfile, append = T) cat("Iterations used = ", data.defaults()$start+globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="", file = globals()$outfile, append = T) cat("Thinning interval =", data.defaults()$thin, "\n", file = globals()$outfile, append = T) cat("Sample size per chain =", globals()$Niters, "\n", file = globals()$outfile, append = T) cat("\nFraction in 1st window =", func.defaults()$frac1, file = globals()$outfile, append = T) cat("\nFraction in 2nd window =", func.defaults()$frac2, "\n\n", file = globals()$outfile, append = T) gcd <- matrix(nrow = globals()$Nchains, ncol = globals()$Nparms+3) for(i in 1:globals()$Nchains) { gcd[i, 1:3] <- c(paste(globals()$filenames[data.defaults()$chains[i]], " "), paste(rep("=", nchar( globals()$filenames[data.defaults()$chains[i]])), collapse=""), "") for(j in 1:(globals()$Nparms)) { gcd[i, j+3] <- geweke.cd(globals()$parmsmat[, j, i]) } } gcd <- t(gcd) gcd[-1:-3,] <- format(signif(as.numeric(gcd[-1:-3,]), digits=func.defaults()$digits)) # # Calculate how many columns will fit on page - if more columns # are needed than will fit, then spread over 2 or more tables # cwidth <- 80 - (max(nchar(globals()$Labels)) + 6) colsize <- apply(nchar(gcd), 2, max) + 5 tmp <- cumsum(colsize) nr <- list() nr[[1]] <- c(1:length(colsize))[tmp <= cwidth] i<-1 while(nr[[i]][length(nr[[i]])] != length(colsize)) { tmp <- tmp + (cwidth*i - tmp[nr[[i]][length(nr[[i]])]]) nr[[i+1]] <- c(1:length(colsize))[tmp > (cwidth*i) & tmp <= (cwidth*(i+1))] i <- i+1 } ntabs <- length(nr) for(n in 1:ntabs) { gcd.table <- pretty.print(matrix(c(globals()$parmnames, pretty.print( gcd[, nr[[n]]], print.it=F, hsep=" ", vsep="",csep=" ")), ncol=2), hsep=" | ", vsep="-", csep="-+-", print.it=F) cat(gcd.table) cat(gcd.table, file=globals()$outfile, append=T) } while(T) { cat("\nGeweke Plots Menu\n") cat( "*****************\n") geweke.menu <- c("Plot Z-scores", "Return to Diagnostics Menu") pick <- menu(geweke.menu) switch(pick, { if(is.numeric(func.defaults()$geweke.bin) & globals()$Niters<50+func.defaults()$geweke.bin) { cat("\n******* Error: *******") cat("\nNot enough iterations available for plot\n") } else { if(is.numeric(func.defaults()$geweke.max) & globals()$Niters<50+func.defaults()$geweke.max) { cat("\n******* Error: *******") cat("\nNot enough iterations available for plot\n") } else { open.graphics() args.list <- list() args.list[[1]] <- numeric() args.list[[2]] <- numeric() args.list[[3]] <- numeric() args.list[[4]] <- numeric() names(args.list) <- c("this.parm", "this.chain", "scale", "title.scale") setup.plots(globals()$Nparms, nplots=1, sepplot=T, which.plot = "geweke.plot", args.list, extra.title = "Geweke's Convergence Diagnostic") } } }, break ) } } # # ----------------------------------------------------------------------- # geweke.plot _ function(this.parm, this.chain, scale, title.scale) { # # geweke.plot -- plots repreated Geweke convergence diagnostic # (Z-score) against iteration # v <- globals()$Labels[data.defaults()$parms] cn <- globals()$filenames[data.defaults()$chains] if(is.numeric(func.defaults()$geweke.bin)) { if(is.numeric(func.defaults()$geweke.max)) { bin.option <- 1 # default option } else { bin.option <- 2 # use bin size } } else { bin.option <- 3 # use max number of bins } # number of bins: nbin <- switch(bin.option, min(floor((globals()$Niters-50)/func.defaults()$geweke.bin), func.defaults()$geweke.max), floor((globals()$Niters-50)/func.defaults()$geweke.bin), func.defaults()$geweke.max ) # bin width: bin <- switch(bin.option, floor((globals()$Niters-50)/nbin), func.defaults()$geweke.bin, floor((globals()$Niters-50)/nbin) ) gcd <- numeric(nbin+1) gcd[1] <- geweke.cd(globals()$parmsmat[, this.parm, this.chain]) if(nbin>2) { for(n in 2:nbin) { gcd[n] <- geweke.cd(globals()$parmsmat[(((n-1)*bin)+1):globals()$Niters, this.parm, this.chain]) } } gcd[nbin+1] <- geweke.cd(globals()$parmsmat[(globals()$Niters-49):globals()$Niters, this.parm, this.chain]) i <- data.defaults()$start+globals()$out[[1]][1, "iter"]-1 # number of first # iteration title.scale<-switch(func.defaults()$ps.plot, title.scale, title.scale*.85, title.scale*.7) scale<-switch(func.defaults()$ps.plot, scale, scale*.8, scale*.6) ymin <- min(min(gcd), -1.96) ymax <- max(max(gcd), 1.96) + 0.4*(max(max(gcd), 1.96)-ymin) lim <- range(pretty(c(ymin, ymax))) plot(c((0:(nbin-1))*bin, globals()$Niters-49)+i, gcd, type = "p", xlab = "First iteration in segment", ylab = "Z-score", cex = scale, ylim=lim, pch=4, axes=F) axis(1, at=pretty(c((1:nbin)*bin, globals()$Niters-49)+i, 3), cex = scale, labels=my.format(pretty(c((1:nbin)*bin, globals()$Niters-49)+i, 3), scientific=c(-func.defaults()$scientific, func.defaults()$scientific))) par(crt=90, srt=90) axis(2, at=pretty(c(ymin, ymax), 2), cex = scale, labels=my.format(pretty(c(ymin, ymax), 2), scientific=c(-func.defaults()$scientific, func.defaults()$scientific))) par(crt=0, srt=0) box() abline(h=1.96, lty=2) abline(h=-1.96, lty=2) if(globals()$Nchains>1) { title(main = paste(v[this.parm], " (", cn[this.chain], ")", sep=""), cex=title.scale) } else { title(main = paste(v[this.parm], sep=""), cex=title.scale) } } # # ----------------------------------------------------------------------- # geweke.cd _ function(x) { # # geweke.cd -- computes Geweke convergence diagnostic # (Modified by NGB 2/4/97 to fix error in specifying window widths) # # x is vector representing one G.S. chain # calls geweke.power a <- length(x) power1 <- geweke.power(x[1:(a*func.defaults()$frac1)]) power2 <- geweke.power(x[(a - (a*func.defaults()$frac2) + 1):a]) gbara <- mean(x[1:(a*func.defaults()$frac1)]) gbarb <- mean(x[(a - (a*func.defaults()$frac2) + 1):a]) # Error fixed 2/4/97 (NGB): changed a/2 to a*func.defaults()$frac1 & a/5 # to a*func.defaults()$frac2 cd <- (gbara - gbarb)/sqrt((power1/(a*func.defaults()$frac1)) + (power2/(a*func.defaults()$frac2))) cd } # # ----------------------------------------------------------------------- # geweke.power _ function(x) { # # geweke.power # a <- length(x) nspans <- sqrt(a)/0.3 + 1 # spans parm for smoothing periodogram if(nspans >= a) { # spans is longer than time series nspans <- a - 1 } pgram1 <- spec.pgram(x, spans = nspans, demean = T, detrend = F, plot = F) power <- (10^(pgram1$spec[1]/10)) #changed 03/20/94 # power <- 2 * pi * (10^(pgram1$spec[1]/10)) # spectral density converted from decibels nse <- sqrt(power/a) power } # # ----------------------------------------------------------------------- # geweke.nse _ function(x) { # # geweke.nse # # Author: Kate Cowles # a <- length(x) nspans <- sqrt(a)/0.3 + 1 # spans parm for smoothing periodogram 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 <- (10^(pgram1$spec[1]/10)) # power <- 2 * pi * (10^(pgram1$spec[1]/10)) # spectral density converted from decibels nse <- sqrt(power/a) nse } # # ----------------------------------------------------------------------- # gandr _ function() { # # gandr -- Calculate Gelman and Rubin convergence diagnostic for CODA system # # Author: Kate Cowles # Modified by: Nicky Best # a <- length(globals()$parmsmat[, 1, 1]) if(globals()$Nchains == 1) { cat("\nError: Gelman and Rubin convergence diagnostic requires more than one chain.\n\n") } else { if(all(globals()$parmsmat[(a/2 + 1):a, 1, 1] == globals()$parmsmat[(a/2 + 1):a, 1, 2])) { cat("\nError: 2nd half of all chains are identical\nGelman and Rubin convergence diagnostic requires 2 or more different chains.\n\n") } else { cat("\nGELMAN AND RUBIN 50% AND 97.5% SHRINK FACTORS:") cat("\n============================================== \n\n") cat("Iterations used for diagnostic = ", ceiling(data.defaults()$start+0.5*(data.defaults()$end-data.defaults()$start)) , ":", data.defaults()$end, "\n", sep="") cat("Thinning interval =", data.defaults()$thin, "\n") cat("Sample size per chain =", globals()$Niters, "\n\n") cat("\nGELMAN AND RUBIN 50% AND 97.5% SHRINK FACTORS:", file = globals()$outfile, append = T) cat("\n============================================== \n\n", file = globals()$outfile, append = T) cat("Iterations used for diagnostic = ", ceiling(data.defaults()$start+0.5*(data.defaults()$end-data.defaults()$start)) , ":", data.defaults()$end, "\n", sep="", file = globals()$outfile, append = T) cat("Thinning interval =", data.defaults()$thin, "\n", file = globals()$outfile, append = T) cat("Sample size per chain =", globals()$Niters, "\n\n", file = globals()$outfile, append = T) shrink <- matrix(c("Point est.", "97.5% quantile", "==========", "==============", "", ""), nrow=2) v <- globals()$Labels[data.defaults()$parms] which.pos <- func.defaults()$pos.parmsmat pos <- rep(FALSE, length(v)) if(length(which.pos) > 1 || (length(which.pos) == 1 && which.pos > 0)) { for(k in 1:length(which.pos)) { pos[which.pos[k]] <- TRUE } } which.prop <- func.defaults()$prop.parmsmat prop <- rep(FALSE, length(v)) if(length(which.prop) > 1 || (length(which.prop) == 1 && which.prop > 0)) { for(k in 1:length(which.prop)) { prop[which.prop[k]] <- TRUE } } for(i in 1:globals()$Nparms) { r <- globals()$parmsmat[,i,] if(pos[i]){ shrink <- cbind(shrink, gpar.log(r)$confshrink) } else { if(prop[i]){ shrink <- cbind(shrink, gpar.logit(r)$confshrink) } else { shrink <- cbind(shrink, gpar(r)$confshrink) } } } shrink <- t(shrink) shrink[-1:-3,] <- format(signif(as.numeric(shrink[-1:-3,]), digits=func.defaults()$digits)) shrink.table <- pretty.print(matrix(c(globals()$parmnames, pretty.print(shrink, print.it=F, hsep=" ", vsep="", csep=" ")), ncol=2), hsep=" | ", vsep="-", csep="-+-", print.it=F) cat(shrink.table) cat(shrink.table, file=globals()$outfile, append=T) while(T) { cat("\nGelman & Rubin Plots Menu\n") cat( "*************************\n") gr.menu <- c("Trace Plots", "Shrink Factor Plots", "Return to Diagnostics Menu") pick <- menu(gr.menu) switch(pick, { open.graphics() args.list <- list() args.list[[1]] <- numeric() args.list[[2]] <- numeric() args.list[[3]] <- numeric() args.list[[4]] <- numeric() args.list[[5]] <- pos args.list[[6]] <- prop names(args.list) <- c("this.parm", "this.chain", "scale", "title.scale", "pos", "prop") setup.plots(globals()$Nparms, nplots=1, sepplot=F, which.plot = "gandr.trace.plot", args.list, extra.title="Traces plus Gelman & Rubin Shrink Factors") }, { if(is.numeric(func.defaults()$gr.bin) & globals()$Niters<=0+func.defaults()$gr.bin) { cat("\n******* Error: *******") cat("\nNot enough iterations available for plot\n") } else { if(is.numeric(func.defaults()$gr.max) & globals()$Niters<50+func.defaults()$gr.max) { cat("\n******* Error: *******") cat("\nNot enough iterations available for plot\n") } else { open.graphics() args.list <- list() args.list[[1]] <- numeric() args.list[[2]] <- numeric() args.list[[3]] <- numeric() args.list[[4]] <- numeric() args.list[[5]] <- pos args.list[[6]] <- prop names(args.list) <- c("this.parm", "this.chain", "scale", "title.scale", "pos", "prop") setup.plots(globals()$Nparms, nplots=1, sepplot=F, which.plot = "gandr.shrink.plot", args.list, extra.title="Gelman & Rubin Shrink Factors") } } }, break ) } } } } # # ----------------------------------------------------------------------- # gandr.trace.plot _ function(this.parm, this.chain, scale, title.scale, pos, prop) { # # gandr.trace.plot -- function for plotting multiple chains on same graph and # printing Gelman & Rubin's convergence diagnostic for each # variable - called by setup.plots() # # Author: Nicky Best # title.scale<-switch(func.defaults()$ps.plot, title.scale*.8, title.scale*.7, title.scale*.65) scale<-switch(func.defaults()$ps.plot, scale, scale*.8, scale*.6) r <- globals()$parmsmat[,this.parm,] lim <- range(pretty(c(min(r), max(r)))) plot(globals()$parmsmat.iter, r[, 1], type = "l", ylim = lim, xlab = "Iteration", ylab = globals()$Labels[data.defaults()$parms[this.parm]], cex = scale, axes=F) axis(1, at=pretty(globals()$parmsmat.iter,3), cex = scale, labels=my.format(pretty(globals()$parmsmat.iter,3), scientific=c(-func.defaults()$scientific, func.defaults()$scientific))) par(crt=90, srt=90) axis(2, at=pretty(c(min(r), max(r)), 2), cex = scale, labels=my.format(pretty(c(min(r), max(r)), 2), scientific=c(-func.defaults()$scientific, func.defaults()$scientific))) par(crt=0, srt=0) box() for(j in 2:globals()$Nchains) { lines(r[, j], lty = j) } if(pos[this.parm]) { shrink <- gpar.log(r)$confshrink } else { if(prop[this.parm]) { shrink <- gpar.logit(r)$confshrink } else { shrink <- gpar(r)$confshrink } } title(main = paste("Median =", format(round(shrink[1], 2)), ",", "97.5% =", format(round(shrink[2], 2))), cex=title.scale) } # # ----------------------------------------------------------------------- # gandr.shrink.plot _ function(this.parm, this.chain, scale, title.scale, pos, prop) { # # gandr.shrink.plot -- function for plotting Gelman & Rubin's shrink factors # against iteration number for each variable - called by # setup.plots() # # Author: Nicky Best # r <- globals()$parmsmat[,this.parm,] # Want 50 iterations in first bin; default then takes bins of size 10 unless # this gives more than 50 bins, in which case, bin size = (Niters-50)/100. if(is.numeric(func.defaults()$gr.bin)) { if(is.numeric(func.defaults()$gr.max)) { bin.option <- 1 # default option } else { bin.option <- 2 # use bin size } } else { bin.option <- 3 # use max number of bins } # number of bins: nbin <- switch(bin.option, min(floor((globals()$Niters-50)/func.defaults()$gr.bin), func.defaults()$gr.max), floor((globals()$Niters-50)/func.defaults()$gr.bin), func.defaults()$gr.max ) # bin width: bin <- switch(bin.option, floor((globals()$Niters-50)/nbin), func.defaults()$gr.bin, floor((globals()$Niters-50)/nbin) ) shrink <- matrix(ncol=2, nrow=nbin+1) if(pos[this.parm]) { shrink[1,] <- gpar.log(r[1:50,])$confshrink if(nbin>2) { for(n in 2:nbin) { shrink[n,] <- gpar.log(r[1:(50+(n-1)*bin),])$confshrink } } shrink[nbin+1,] <- gpar.log(r)$confshrink } else { if(prop[this.parm]) { shrink[1,] <- gpar.logit(r[1:50,])$confshrink if(nbin>2) { for(n in 2:nbin) { shrink[n,] <- gpar.logit(r[1:(50+(n-1)*bin),])$confshrink } } shrink[nbin+1,] <- gpar.logit(r)$confshrink } else { shrink[1,] <- gpar(r[1:50,])$confshrink if(nbin>2) { for(n in 2:nbin) { shrink[n,] <- gpar(r[1:(50+(n-1)*bin),])$confshrink } } shrink[nbin+1,] <- gpar(r)$confshrink } } i <- data.defaults()$start+globals()$out[[1]][1, "iter"]-2 # number of first # iteration -1 title.scale<-switch(func.defaults()$ps.plot, title.scale, title.scale*.85, title.scale*.7) scale<-switch(func.defaults()$ps.plot, scale, scale*.8, scale*.6) if(length(shrink[!is.na(shrink)])==0) { cat("\n******* Error: *******\nCannot compute Gelman & Rubin's diagnostic for any chain segments\nThis indicates convergence failure ==> Run chains for more iterations\n") } else { ymin <- min(min(shrink[!is.na(shrink)]), 1) ymax <- max(max(shrink[!is.na(shrink)]), 1) + 0.4*(max(max(shrink[!is.na(shrink)]), 1)-ymin) lim <- range(pretty(c(ymin, ymax))) plot(c(50, 50+(1:nbin)*bin)+i, shrink[,1], type = "l", xlab = "Last iteration in segment", ylab = "Shrink factor", cex = scale, ylim=lim, xlim=c(50+i, 50+nbin*bin+i), axes=F) axis(1, at=pretty(c(50, 50+(1:nbin)*bin)+i, 3), cex = scale, labels=my.format(pretty(c(50, 50+(1:nbin)*bin)+i, 3), scientific=c(-func.defaults()$scientific, func.defaults()$scientific))) par(crt=90, srt=90) axis(2, at=pretty(c(ymin, ymax), 2), cex = scale, labels=my.format(pretty(c(ymin, ymax), 2), scientific=c(-func.defaults()$scientific, func.defaults()$scientific))) par(crt=0, srt=0) box() lines(c(50, 50+(1:nbin)*bin)+i, shrink[,2], lty=2) abline(h=1, lty=4) legend(i+50+0.7*((50+nbin*bin+i) - (i+50)), lim[2]+0.02*(lim[2]-lim[1]), c("median", "", "97.5%"), lty=c(1, -1, 2), cex=scale*0.7, bty="n") title(main = globals()$Labels[data.defaults()$parms[this.parm]], cex=title.scale) if(length(shrink[is.na(shrink)])>0) { cat("\n******* Warning: *******\nCould not compute Gelman & Rubin's diagnostic for", length(shrink[is.na(shrink)]), "chain segments\n") } } } # # ----------------------------------------------------------------------- # gpar _ function(r) { # # gpar -- Gelman and Rubin's code # (Modified by NGB 28/8/97 to fix error in correction factor) # alpha <- 0.05 # 95% intervals m <- ncol(r) x <- r[(nrow(r)/2 + 1):nrow(r), ] # 2nd 1/2 of simulated sequences n <- nrow(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 xdot <- as.vector(col.means(x)) s2 <- as.vector(col.vars(x)) W <- mean(s2) B <- n * var(xdot) muhat <- mean(xdot) varW <- var(s2)/m varB <- (B^2 * 2)/(m - 1) covWB <- (n/m) * (cov(s2, xdot^2) - 2 * muhat * cov(s2, xdot)) sig2hat <- ((n - 1) * W + B)/n quantiles <- quantile(as.vector(x), probs = c(0.025, 0.25, 0.5, 0.75, 0.975)) if(W > 1e-08) { # non-degenerate case # 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/(m * n) varpostvar <- (((n - 1)^2) * varW + (1 + 1/m)^2 * varB + 2 * (n - 1) * (1 + 1/m) * covWB)/n^2 post.df <- chisqdf(postvar, varpostvar) post.range <- muhat + sqrt(postvar) * qt(1 - alpha/2, post.df) * c(-1, 0, 1) # Estimated potential scale reduction (that would be achieved by # continuing simulations forever) has two components: an estimate and # an approx. 97.5% 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. # NOTE: original correction factor used by Gelman & Rubin (1992) and # and used in the original CODA (version 0.3) code was sqrt(df/(df-2)). # This is incorrect and can lead to estimates of the scale reduction # factor which are infinite or negative. The correct correction factor # is sqrt((df+3)/(df+1)). See Brooks S.P. & Gelman A. (1998) "General # Methods for Monitoring Convergence of Iterative Simulations", # Journal of Computational and Graphical Statistics (to appear) - also # available from the MCMC preprint service # (http: //www.stats.bris.ac.uk/MCMC) # # 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 <- chisqdf(W, varW) confshrink.range <- sqrt((c(postvar/W, (n - 1)/n + (1 + 1/m) * ( 1/n) * (B/W) * qf(0.975, m - 1, varlo.df)) * (post.df+3))/(post.df +1)) # correction factor edited by NGB 28/8/97 list(post = post.range, quantiles = quantiles, confshrink = confshrink.range) } else { # degenerate case: all entries in "data matrix" are identical list(post = muhat * c(1, 1, 1), quantiles = quantiles, confshrink = c(1, 1)) } } # # ----------------------------------------------------------------------- # gpar.log _ function(r) { gp <- gpar(log(r)) list(post = exp(gp$post), quantiles = exp(gp$quantiles), confshrink = gp$confshrink) } # # ----------------------------------------------------------------------- # gpar.logit _ function(r) { gp <- gpar(logit(r)) list(post = invlogit(gp$post), quantiles = invlogit(gp$quantiles), confshrink = gp$confshrink) } # # ----------------------------------------------------------------------- # logit _ function(x) { log(x/(1 - x)) } # # ----------------------------------------------------------------------- # invlogit _ function(x) { 1/(1 + exp(.Uminus(x))) } # # ----------------------------------------------------------------------- # col.means _ function(mat) { # # Called by gpar() for calculating Gelman & Rubin diagnostic: # ones <- matrix(1, nrow = 1, ncol = nrow(mat)) ones %*% mat/nrow(mat) } # # ----------------------------------------------------------------------- # col.vars _ function(mat) { # # Called by gpar() for calculating Gelman & Rubin diagnostic: # means <- col.means(mat) col.means(mat * mat) - means * means } # # ----------------------------------------------------------------------- # cov _ function(a, b) { # # Called by gpar() for calculating Gelman & Rubin diagnostic: # m <- length(a) ((mean((a - mean(a)) * (b - mean(b)))) * m)/(m - 1) } # # ----------------------------------------------------------------------- # chisqdf _ function(A, varA) { # # Called by gpar() for calculating Gelman & Rubin diagnostic: # 2 * (A^2/varA) } # # ----------------------------------------------------------------------- # autocorr _ function() { # # autocorr -- Calculate Autocorrelations for CODA system # # Author: Kate Cowles # Modified by: Nicky Best # cat("\nLAGS AND AUTOCORRELATIONS WITHIN EACH CHAIN:") cat("\n============================================\n\n") cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="") cat("Thinning interval =", data.defaults()$thin, "\n") cat("Sample size per chain =", globals()$Niters, "\n\n") cat("\nLAGS AND AUTOCORRELATIONS WITHIN EACH CHAIN:", file = globals()$outfile, append = T) cat("\n============================================\n\n", file = globals()$outfile, append = T) cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="", file = globals()$outfile, append = T) cat("Thinning interval =", data.defaults()$thin, "\n", file = globals()$outfile, append = T) cat("Sample size per chain =", globals()$Niters, "\n\n", file = globals()$outfile, append = T) ac.mat <- matrix(ncol=4, nrow=2) for(i in 1:globals()$Nchains) { for(j in 1:globals()$Nparms) { temp <- acf(globals()$parmsmat[, j, i], lag.max = 50, plot = F) ac.mat <- rbind(ac.mat, temp$acf[c(2, 6, 11, 51)]) } } ac.mat[-1:-2,] <- format(signif(as.numeric(ac.mat[-1:-2,]), digits=func.defaults()$digits)) if(nrow(ac.mat) == 3) { lab.width <- nchar(ac.mat[-1:-2,]) } else { lab.width <- apply(nchar(ac.mat[-1:-2,]), 2, max) } extra.space <- lab.width - c(5, 5, 6, 6) extra.space[extra.space < 0] <- 0 space <- numeric() for(i in 1:4) { space[i] <- paste(rep(" ", extra.space[i]), collapse="") } ac.mat[1:2, ] <- matrix(c(paste(paste("Lag", c(1, 5, 10, 50)), space, sep=""), "=====", "=====", "======", "======"), ncol=4, byrow=T) ac.chains <- matrix(c("Chain", "====="), ncol=1) ac.vars <- matrix(c("Variable", "========", globals()$Labels[data.defaults()$parms]), ncol=1) # # Calculate how many columns will fit on page - if more columns # are needed than will fit, then spread over 2 or more tables # cwidth <- 80 - (max(nchar(globals()$Labels)) + 6) colsize <- apply(nchar(ac.mat), 2, max) + 5 tmp <- cumsum(colsize) nr <- list() nr[[1]] <- c(1:length(colsize))[tmp <= cwidth] i<-1 while(nr[[i]][length(nr[[i]])] != length(colsize)) { tmp <- tmp + (cwidth*i - tmp[nr[[i]][length(nr[[i]])]]) nr[[i+1]] <- c(1:length(colsize))[tmp > (cwidth*i) & tmp <= (cwidth*(i+1))] i <- i+1 } ntabs <- length(nr) for(n in 1:ntabs) { ac.table <- cbind(pretty.print(ac.chains, print.it=F, hsep="", vsep="", csep=""), pretty.print(ac.vars[1:2, 1], print.it=F, hsep=" ", vsep="", csep=" "), pretty.print(ac.mat[1:2, nr[[n]]], print.it=F, hsep=" ", vsep="", csep=" ", min.colwidth = 9)) for(j in 1:globals()$Nchains) { ac.table <- rbind(ac.table, matrix(c(globals()$filenames[data.defaults()$chains[j]], pretty.print(ac.vars[-1:-2,1], print.it=F, hsep=" ", vsep="", csep=" "), pretty.print(matrix(ac.mat[((j-1)*globals()$Nparms+3): (j*globals()$Nparms+2), nr[[n]]], ncol= length(nr[[n]])), print.it=F, hsep=" ", vsep="", csep=" ", min.colwidth = 9)), nrow=1)) } ac.table <- pretty.print(ac.table, hsep=" | ", vsep="-", csep="-+-", print.it=F) cat(ac.table) cat(ac.table, file= globals()$outfile, append = T) } while(T) { cat("\nAutocorrelation Plots Menu\n") cat( "**************************\n") autocor.menu <- c("Plot autocorrelations", "Return to Diagnostics Menu") pick <- menu(autocor.menu) switch(pick, { open.graphics() args.list <- list() args.list[[1]] <- numeric() args.list[[2]] <- numeric() args.list[[3]] <- numeric() args.list[[4]] <- numeric() names(args.list) <- c("this.parm", "this.chain", "scale", "title.scale") setup.plots(globals()$Nparms, nplots=1, sepplot=T, which.plot = "autocorr.plot", args.list, extra.title = "Autocorrelations") }, break ) } } # # ----------------------------------------------------------------------- # autocorr.plot _ function(this.parm, this.chain, scale, title.scale) { # # autocorr.plot -- function to plot autocorrelations within each chain for # each parameter - called by setup.plots() # # Author: Nicky Best # title.scale<-switch(func.defaults()$ps.plot, title.scale, title.scale*.85, title.scale*.7) scale<-switch(func.defaults()$ps.plot, scale, scale*.8, scale*.6) temp <- acf(globals()$parmsmat[, this.parm, this.chain], lag.max = 50, plot = F) plot(temp$acf, type = "h", ylab = "Autocorrelation", xlab = "Lag", cex = scale, ylim=c(-1, 1)) title(paste(globals()$Labels[data.defaults()$parms[this.parm]], ": ", globals()$filenames[data.defaults()$chains[this.chain]], sep=""), cex=title.scale) } # # ----------------------------------------------------------------------- # parmcorr _ function() { # # parmcorr -- Cross-correlations for CODA system # # Author: Kate Cowles # Modified by: Nicky Best # cat("\nCROSS-CORRELATION MATRIX:") cat("\n=========================\n\n") cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="") cat("Thinning interval =", data.defaults()$thin, "\n") cat("Sample size per chain =", globals()$Niters, "\n\n") cat("\nCROSS-CORRELATION MATRIX:", file = globals()$outfile, append = T) cat("\n=========================\n\n", file = globals()$outfile, append = T) cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="", file = globals()$outfile, append = T) cat("Thinning interval =", data.defaults()$thin, "\n", file = globals()$outfile, append = T) cat("Sample size per chain =", globals()$Niters, "\n\n", file = globals()$outfile, append = T) if(globals()$Nchains>1) { if(func.defaults()$combine.corr) { cat("Pooling over chains: ") cat("Pooling over chains: ", file=globals()$outfile, append = T) cat(globals()$filenames[data.defaults()$chains[1]], "\n") cat(globals()$filenames[data.defaults()$chains[1]], "\n", file=globals()$outfile, append = T) for(i in 2:globals()$Nchains) { cat(" ", globals()$filenames[data.defaults()$chains[i]], "\n", sep="") cat(" ", globals()$filenames[data.defaults()$chains[i]], "\n", sep="", file=globals()$outfile, append = T) } cat("\n") cat("\n", file=globals()$outfile, append = T) } } if(func.defaults()$combine.corr) { Nchains.corr <- 1 } else Nchains.corr <- globals()$Nchains for(j in 1:Nchains.corr) { if(func.defaults()$combine.corr) { combine.parmsmat <- globals()$parmsmat[, ,1] for(k in 2:globals()$Nchains) { combine.parmsmat <- rbind(combine.parmsmat, globals()$parmsmat[, ,k]) } corr.mat <- format(signif(cor(combine.parmsmat), digits=func.defaults()$digits)) } else { corr.mat <- format(signif(cor(globals()$parmsmat[, ,j]), digits=func.defaults()$digits)) } dim(corr.mat) <- c(globals()$Nparms, globals()$Nparms) if(globals()$Nparms > 1) { for(rw in 1:(globals()$Nparms-1)) { for(cl in (rw+1):globals()$Nparms) { corr.mat[rw, cl] <- "" } } } corr.mat <- rbind(matrix(c(globals()$Labels[data.defaults()$parms], rep("", 2*globals()$Nparms)), byrow=T, nrow=3), corr.mat) if(!func.defaults()$combine.corr) { cat("\nChain: ", globals()$filenames[data.defaults()$chains[j]], "\n", paste(c("=======", rep("=", nchar(globals()$filenames[data.defaults()$chains[j]]))), collapse=""), "\n\n", sep="") cat("\nChain: ", globals()$filenames[data.defaults()$chains[j]], "\n", paste(c("=======", rep("=", nchar(globals()$filenames[data.defaults()$chains[j]]))), collapse=""), "\n\n", sep="", file=globals()$outfile, append=T) } corr.table <- list() # # Calculate how many columns will fit on page - if more columns # are needed than will fit, then spread over 2 or more tables # cwidth <- 80 - (max(nchar(globals()$Labels)) + 6) colsize <- apply(nchar(corr.mat), 2, max) + 5 tmp <- cumsum(colsize) nr <- list() nr[[1]] <- c(1:length(colsize))[tmp <= cwidth] i<-1 while(nr[[i]][length(nr[[i]])] != length(colsize)) { tmp <- tmp + (cwidth*i - tmp[nr[[i]][length(nr[[i]])]]) nr[[i+1]] <- c(1:length(colsize))[tmp > (cwidth*i) & tmp <= (cwidth*(i+1))] i <- i+1 } ntabs <- length(nr) for(n in 1:ntabs) { corr.table[[n]] <- pretty.print(matrix(c(globals()$parmnames, pretty.print(corr.mat[, nr[[n]]], print.it=F, hsep=" ", vsep="",csep=" ")), ncol=2), hsep=" | ", vsep="-", csep="-+-", print.it=F) cat(corr.table[[n]]) cat(corr.table[[n]], file=globals()$outfile, append=T) } } while(T) { cat("\nCross Correlation Plots Menu\n") cat( "****************************\n") autocor.menu <- c("Plot Cross Correlations", "Return to Diagnostics Menu") pick <- menu(autocor.menu) switch(pick, { open.graphics() corr.list <- list() for(j in 1:Nchains.corr) { if(func.defaults()$combine.corr) { corr.list[[j]] <- as.matrix(cor(combine.parmsmat)) } else { corr.list[[j]] <- as.matrix(cor(globals()$parmsmat[, , j])) } } args.list <- list() args.list[[1]] <- corr.list args.list[[2]] <- numeric() args.list[[3]] <- 1 args.list[[4]] <- numeric() args.list[[5]] <- numeric() args.list[[6]] <- numeric() names(args.list) <- c("corr.list", "this.parm", "this.chain", "scale", "title.scale", "mrows") setup.plots(Nparms.tmp=1, nplots=1, sepplot=T, which.plot = "parmcorr.plot", args.list, one.page = F, extra.title= "Cross Correlations") }, break ) } } # # ----------------------------------------------------------------------- # parmcorr.plot _ function(corr.list, this.chain, this.parm, scale, title.scale, mrows) { # # parmcorr.plot -- function to plot autocorrelations within each chain for # each parameter - called by setup.plots() # # Author: Nicky Best # title.scale<-switch(func.defaults()$ps.plot, title.scale, title.scale*.85, title.scale*.7) scale<-switch(func.defaults()$ps.plot, scale, scale*.8, scale*.6) scale1<-switch(func.defaults()$ps.plot, scale*.6, scale*.7, scale*.7) par(xaxt="n", yaxt="n") plot(0, 0, type="n", xlim=c(0, globals()$Nparms), ylim=c(0, globals()$Nparms), xlab="", ylab="", cex=scale) par(adj=1) axis(2, at=c(globals()$Nparms:1)-0.5, labels=globals()$short.labels[data.defaults()$parms], cex=scale) par(crt=90, srt=90) axis(1, at=c(1:globals()$Nparms)-0.5, labels=globals()$short.labels[data.defaults()$parms], cex=scale) par(xaxt="s", yaxt="s", crt=0, srt=0, adj=0.5) for(cl in 1:globals()$Nparms) { for(rw in 1:(globals()$Nparms-cl+1)) { dens <- switch(abs(corr.list[[this.chain]][globals()$Nparms-rw+1, cl]*10) %/% 2 + 1, 0, 10, 25, 40, 65, 100 ) # # Fudge to ensure density on diagonal = 1.0, # as rounding errors may occur in switch: # if(cl == globals()$Nparms-rw+1) dens <- 100 # ang <- ifelse(corr.list[[this.chain]][globals()$Nparms-rw+1, cl] < 0, 315, 45) polygon(c(cl-1, cl-1, cl, cl, cl-1), c(rw-1, rw, rw, rw-1, rw-1), density = dens, angle = ang) } } leg.scale <- c(0.35, 0.675) if(mrows ==1) { legend(c(leg.scale[1] * globals()$Nparms, leg.scale[2] * globals()$Nparms), c(0.99 * globals()$Nparms, 0.75 * globals()$Nparms), legend = c("|r| = [0.0, 0.2)", "", "|r| = [0.2, 0.4)", "", "|r| = [0.4, 0.6)"), density = c(0.1, -99, 10, -99, 25), angle = rep(45, 5), cex = scale1*0.85, bty = "n") legend(c(leg.scale[2] * globals()$Nparms, 1.02 * globals()$Nparms), c(0.99 * globals()$Nparms, 0.75 * globals()$Nparms), legend = c( "|r| = [0.6, 0.8)", "", "|r| = [0.8, 1.0)", "", "|r| = 1.0"), density = c(40, -99, 65, -99, 100), angle = rep(45, 5), cex = scale1*0.85, bty = "n") par(adj = 0) text(0.5 * (leg.scale[1] + leg.scale[2]) * globals()$Nparms, 1.025 * globals()$Nparms, "Gradient=sign of correlation", cex = scale1*0.85) } else { legend(c(leg.scale[1] * globals()$Nparms, leg.scale[2] * globals()$Nparms), c(0.97 * globals()$Nparms, 0.8 * globals()$Nparms), legend = c("|r| = [0.0, 0.2)", "|r| = [0.2, 0.4)", "|r| = [0.4, 0.6)"), density = c(0.1, 10, 25), angle = rep(45, 3), cex = scale1*0.8, bty = "n") legend(c(leg.scale[2] * globals()$Nparms, 1.02 * globals()$Nparms), c(0.97 * globals()$Nparms, 0.8 * globals()$Nparms), legend = c( "|r| = [0.6, 0.8)", "|r| = [0.8, 1.0)", "|r| = 1.0"), density = c(40, 65, 100), angle = rep(45, 3), cex = scale1*0.8, bty = "n") par(adj = 0) text(0.5 * (leg.scale[1] + leg.scale[2]) * globals()$Nparms, globals()$Nparms, "Gradient=sign of correlation", cex = scale1*0.8) } par(adj=0.5) if(globals()$Nchains>1) { title(globals()$filenames[data.defaults()$chains[this.chain]], cex=title.scale) } } # # ----------------------------------------------------------------------- # heidel _ function() { # # Heidelberger and Welch diagnostic # # Author: Kate Cowles # Modified by: Nicky Best # eps <- func.defaults()$halfwidth cat("\nHEIDELBERGER AND WELCH STATIONARITY AND INTERVAL HALFWIDTH TESTS:") cat("\n=================================================================\n\n") cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="") cat("Thinning interval =", data.defaults()$thin, "\n") cat("Sample size per chain =", globals()$Niters, "\n\n") cat("Precision of halfwidth test =", func.defaults()$halfwidth, "\n\n") cat("\nHEIDELBERGER AND WELCH STATIONARITY AND INTERVAL HALFWIDTH TESTS:", file = globals()$outfile, append = T) cat("\n=================================================================\n\n", file = globals()$outfile, append = T) cat("Iterations used = ", data.defaults()$start +globals()$out[[1]][1, "iter"]-1, ":", data.defaults()$end+globals()$out[[1]][1, "iter"]-1, "\n", sep="", file = globals()$outfile, append = T) cat("Thinning interval =", data.defaults()$thin, "\n", file = globals()$outfile, append = T) cat("Sample size per chain =", globals()$Niters, "\n\n", file = globals()$outfile, append = T) cat("Precision of halfwidth test =", func.defaults()$halfwidth, "\n\n", file = globals()$outfile, append = T) HW.mat <- NULL for(i in 1:globals()$Nchains) { HW.mat <- matrix(c("Stationarity ", " test", "============", "", "# of iters. ", " to keep", "===========", "", "# of iters. ", "to discard", "===========", "", "C-vonM", " stat.", "======", "", "Halfwidth ", " test", "=========", "", "", "Mean ", "====", "", "", "Halfwidth", "=========", ""), nrow=4) for(j in 1:globals()$Nparms) { Y <- globals()$parmsmat[, j, i] n1 <- length(Y) n <- length(Y) S0 <- geweke.power(Y[(n/2 + 1):n]) passed <- F passed2 <- 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) if(!is.na(I) & I < 0.46) { passed <- T } else { Y <- Y[(n1/10 + 1):n] n <- length(Y) } } S0ci <- geweke.power(Y) halfwidth <- 1.96 * sqrt(S0ci/n) if(!is.na(halfwidth) & abs(halfwidth/ybar) <= eps) { passed2 <- T } if(is.na(I) | is.na(halfwidth) | !passed) { n <- NA nd <- NA passed2 <- NA ybar <- NA halfwidth <- NA } else { nd <- length(globals()$parmsmat[, j, i]) - n passed2 <- ifelse(passed2, "passed", "failed") } passed <- ifelse(passed, "passed", "failed") HW.mat <- rbind(HW.mat, c(passed, n, nd, I, passed2, ybar, halfwidth)) } HW.mat[-1:-4, 4] <- format(signif(as.numeric(HW.mat[-1:-4, 4]), digits=3)) HW.mat[-1:-4, 6] <- format(signif(as.numeric(HW.mat[-1:-4, 6]), digits=func.defaults()$digits)) HW.mat[-1:-4, 7] <- format(signif(as.numeric(HW.mat[-1:-4, 7]), digits=func.defaults()$digits)) HW.table1 <- pretty.print(matrix(c(globals()$parmnames2, pretty.print(HW.mat[,1:4], print.it=F, hsep=" ", vsep="", csep=" ")), ncol=2), hsep=" | ", vsep="-", csep="-+-", print.it=F) HW.table2 <- pretty.print(matrix(c(globals()$parmnames2, pretty.print(HW.mat[,5:7], print.it=F, hsep=" ", vsep="", csep=" ")), ncol=2), hsep=" | ", vsep="-", csep="-+-", print.it=F) cat("\nChain: ", globals()$filenames[data.defaults()$chains[i]], "\n", paste(c("=======", rep("=", nchar(globals()$filenames[data.defaults()$chains[i]]))), collapse=""), "\n\n", sep="") cat("\nChain: ", globals()$filenames[data.defaults()$chains[i]], "\n", paste(c("=======", rep("=", nchar(globals()$filenames[data.defaults()$chains[i]]))), collapse=""), "\n\n", sep="", file=globals()$outfile, append=T) cat(HW.table1) cat(HW.table2) cat(HW.table1, file=globals()$outfile, append = T) cat(HW.table2, file=globals()$outfile, append = T) } } # # ----------------------------------------------------------------------- # display.defaults _ function(data = T, funcs = T) { # # display.defaults -- displays current values of user-controlled parameters # in tabular format # Author: Nicky Best # # if(data) { cat("\nCurrent default settings:") cat("\n=========================\n\n") # # This bit calculates how many variable names will fit onto each row # (maximum space availble = 80 - 26 (for labels & spacing) = 54 columns) # Labels1 <- globals()$Labels[data.defaults()$parms] Labels1[-globals()$Nparms] <- paste(Labels1[-globals()$Nparms], ", ", sep="") tmp <- cumsum(nchar(Labels1)) lr <- list() lr[[1]] <- c(1:globals()$Nparms)[tmp <= 54] rr <- Labels1[lr[[1]]] i<-1 while(rr[length(rr)] != Labels1[globals()$Nparms]) { tmp <- tmp + (54 - sum(nchar(rr))) lr[[i+1]] <- c(1:globals()$Nparms)[tmp > (54*i) & tmp <= (54*(i+1))] rr <- Labels1[lr[[i+1]]] i <- i+1 } nlr <- length(lr) # # This bit calculates how many chain names will fit onto each row # Chains1 <- globals()$filenames[data.defaults()$chains] Chains1[-globals()$Nchains] <- paste(Chains1[-globals()$Nchains], ", ", sep="") tmp <- cumsum(nchar(Chains1)) cr <- list() cr[[1]] <- c(1:globals()$Nchains)[tmp <= 54] rr <- Chains1[cr[[1]]] i<-1 while(rr[length(rr)] != Chains1[globals()$Nchains]) { tmp <- tmp + (54 - sum(nchar(rr))) cr[[i+1]] <- c(1:globals()$Nchains)[tmp > (54*i) & tmp <= (54*(i+1))] rr <- Chains1[cr[[i+1]]] i <- i+1 } ncr <- length(cr) cat("+------------------------------------------------------------------------------+\n") cat("| DATA DEFAULTS |\n") cat("+---------------------+--------------------------------------------------------+\n") cat("| | |\n") cat("| Variables selected: | ") cat(Labels1[lr[[1]]], paste(rep(" ", 54-sum(nchar(Labels1[lr[[1]]]))), collapse=""), " |\n", sep="") if(nlr > 1) { for(i in 2:nlr) { cat("| | ") cat(Labels1[lr[[i]]], paste(rep(" ", 54-sum(nchar(Labels1[lr[[i]]]))), collapse=""), " |\n", sep="") } } cat("| | |\n") cat("| Chains selected: | ") cat(Chains1[cr[[1]]], paste(rep(" ", 54-sum(nchar(Chains1[lr[[1]]]))), collapse=""), " |\n", sep="") if(ncr > 1) { for(i in 2:ncr) { cat(Chains1[cr[[i]]], paste(rep(" ", 54-sum(nchar(Chains1[lr[[i]]]))), collapse=""), " |\n", sep="") } } cat("| | |\n") cat("| Iterations - start: | ") cat(data.defaults()$start+globals()$out[[1]][1, "iter"]-1, paste(rep(" ", 54-nchar(data.defaults()$start)), collapse=""), " |\n", sep="") cat("| end: | ") cat(data.defaults()$end+globals()$out[[1]][1, "iter"]-1, paste(rep(" ", 54-nchar(data.defaults()$end)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Thinning interval: | ") cat(data.defaults()$thin, paste(rep(" ", 54-nchar(data.defaults()$thin)), collapse=""), " |\n", sep="") cat("| | |\n") cat("+---------------------+--------------------------------------------------------+\n") } if(funcs) { # # This bit formats the bandwidth function for printing # funcname <- as.character(func.defaults()$bandwidth)[ length(func.defaults()$bandwidth)] funcn <- nchar(funcname) funcname <- substring(funcname, first = c(1:funcn), last = c(1:funcn)) funcname[funcname == "\n" | funcname == "\t"] <- "" funcname <- paste(funcname[1:funcn], collapse = "") funcn <- nchar(funcname) switch(funcn %/% 54 + 1, funcname <- substring(funcname, first = 1, last = funcn), funcname <- substring(funcname, first = c(1, ( funcn/2 + 1) %/% 1), last = c((funcn/2) %/% 1, funcn)), funcname <- substring(funcname, first = c(1, ( funcn/3 + 1) %/% 1, ((2 * funcn)/3 + 1) %/% 1 ), last = c((funcn/3) %/% 1, ((2 * funcn)/3) %/% 1, funcn)), break) if(funcn %/% 54 + 1 > 3) { # # Truncate function for printing # funcname <- substring(funcname, first = c(1, 55, 109, 164,), last = c(54, 108, 163, 175)) funcname[4] <- paste(funcname[4], "...", sep = "") } # # This bit calculates how many quantiles will fit onto one row # Q1 <- paste(func.defaults()$quantiles * 100, "%", sep="") Q1[-length(Q1)] <- paste(Q1[-length(Q1)], ", ", sep="") tmp <- cumsum(nchar(Q1)) qr <- list() qr[[1]] <- c(1:length(Q1))[tmp <= 54] rr <- Q1[qr[[1]]] i<-1 while(rr[length(rr)] != Q1[length(Q1)]) { tmp <- tmp + (54 - sum(nchar(rr))) qr[[i+1]] <- c(1:length(Q1))[tmp > (54*i) & tmp <= (54*(i+1))] rr <- Q1[qr[[i+1]]] i <- i+1 } nqr <- length(qr) cat("| SUMMARY STATISTIC DEFAULTS |\n") cat("+---------------------+--------------------------------------------------------+\n") cat("| | |\n") cat("| Significant digits: | ") cat(func.defaults()$digits, paste(rep(" ", 54-nchar(func.defaults()$digits)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Combine chains: | ") ans <- ifelse(func.defaults()$combine.stats, "Yes", "No") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Batch size: | ") cat(func.defaults()$batch.size, paste(rep(" ", 54-nchar(func.defaults()$batch.size)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Quantiles: | ") cat(Q1[qr[[1]]], paste(rep(" ", 54-sum(nchar(Q1[qr[[1]]]))), collapse=""), " |\n", sep="") if(nqr > 1) { for(i in 2:nqr) { cat("| | ") cat(Q1[qr[[i]]], paste(rep(" ", 54-sum(nchar(Q1[qr[[i]]]))), collapse=""), " |\n", sep="") } } cat("| | |\n") cat("+---------------------+--------------------------------------------------------+\n") cat("| PLOTTING DEFAULTS |\n") cat("+---------------------+--------------------------------------------------------+\n") cat("| | |\n") cat("| Trace: | ") ans <- ifelse(func.defaults()$trace, "Yes", "No") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Smooth lines: | ") ans <- ifelse(func.defaults()$lowess, "Yes", "No") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Density: | ") ans <- ifelse(func.defaults()$densplot, "Yes", "No") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Bandwidth: | ") cat(funcname[1], paste(rep(" ", 54-nchar(funcname[1])), collapse=""), " |\n", sep="") if(length(funcname) > 1) { for(i in 2:length(funcname)) { cat("| | ") cat(funcname[i], paste(rep(" ", 54-nchar(funcname[i])), collapse=""), " |\n", sep="") } } cat("| | |\n") cat("| Separate plot/chain:| ") ans <- ifelse(func.defaults()$sepplot, "Yes", "No") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Separate page/var.: | ") ans <- ifelse(func.defaults()$onepage, "Yes", "No") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Postscript options: | ") ans <- switch(func.defaults()$ps.plot, "Full page; Landscape", "Full page; Portrait", "Half page; Portrait") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Exponent size: | ") cat(func.defaults()$scientific, paste(rep(" ", 54-nchar(func.defaults()$scientific)), collapse=""), " |\n", sep="") cat("| | |\n") cat("+---------------------+--------------------------------------------------------+\n") cat("| DIAGNOSTICS DEFAULTS |\n") cat("+---------------------+--------------------------------------------------------+\n") cat("| | |\n") cat("| Geweke | |\n") cat("| ------ | |\n") cat("| Window 1 fraction: | ") cat(func.defaults()$frac1, paste(rep(" ", 54-nchar(func.defaults()$frac1)), collapse=""), " |\n", sep="") cat("| Window 2 fraction: | ") cat(func.defaults()$frac2, paste(rep(" ", 54-nchar(func.defaults()$frac2)), collapse=""), " |\n", sep="") cat("| Bin width: | ") if(is.numeric(func.defaults()$geweke.bin)) { cat(func.defaults()$geweke.bin, paste(rep(" ", 54-nchar(func.defaults()$geweke.bin)), collapse=""), " |\n", sep="") } else { cat("N/A", paste(rep(" ", 51), collapse=""), " |\n", sep="") } cat("| Max number of bins: | ") if(is.numeric(func.defaults()$geweke.max)) { cat(func.defaults()$geweke.max, paste(rep(" ", 54-nchar(func.defaults()$geweke.max)), collapse=""), " |\n", sep="") } else { cat("N/A", paste(rep(" ", 51), collapse=""), " |\n", sep="") } cat("| | |\n") cat("| Gelman & Rubin | |\n") cat("| -------------- | |\n") cat("| Bin width: | ") if(is.numeric(func.defaults()$geweke.bin)) { cat(func.defaults()$gr.bin, paste(rep(" ", 54-nchar(func.defaults()$gr.bin)), collapse=""), " |\n", sep="") } else { cat("N/A", paste(rep(" ", 51), collapse=""), " |\n", sep="") } cat("| Max number of bins: | ") if(is.numeric(func.defaults()$gr.max)) { cat(func.defaults()$gr.max, paste(rep(" ", 54-nchar(func.defaults()$gr.max)), collapse=""), " |\n", sep="") } else { cat("N/A", paste(rep(" ", 51), collapse=""), " |\n", sep="") } cat("| | |\n") cat("| Raftery & Lewis | |\n") cat("| --------------- | |\n") cat("| Quantile (q): | ") cat(func.defaults()$q, paste(rep(" ", 54-nchar(func.defaults()$q)), collapse=""), " |\n", sep="") cat("| Precision (+/- r): | ") cat(func.defaults()$r, paste(rep(" ", 54-nchar(func.defaults()$r)), collapse=""), " |\n", sep="") cat("| Probability (s): | ") cat(func.defaults()$s, paste(rep(" ", 54-nchar(func.defaults()$s)), collapse=""), " |\n", sep="") cat("| | |\n") cat("| Cross-correlations | |\n") cat("| ------------------ | |\n") cat("| Combine chains: | ") ans <- ifelse(func.defaults()$combine.corr, "Yes", "No") cat(ans, paste(rep(" ", 54-nchar(ans)), collapse=""), " |\n", sep="") cat("| | |\n") cat("+---------------------+--------------------------------------------------------+\n") } } # # ----------------------------------------------------------------------- # pretty.print _ function(x, ..., hsep = "|", vsep = "-", csep = "+", print.it = T, prefix.width, min.colwidth = .Options$digits, abbreviate.dimnames = T, page.width = .Options$width) { # # pretty.print -- print formatted table # # Author: S-Plus version 3.2 function # Modified by: Nicky Best # # This function is virtually the same as the "print.char.matrix()" # function available in S-Plus version 3.2. It prints a matrix of character # strings. The strings typically have \n's in them to indicate # line breaks and this function will break those strings into separate # lines and put all the lines for each string into one matrix cell. # The matrix cells are separated by the horizontal and vertical separator # characters, hsep and vsep (csep is used where they cross). # This function does not attempt to break up a large matrix into pieces # which fit on a page. # Currently, each call to cat() prints exactly one line. One could # replace the calls to cat with calls to paste(sep="",...) to collect # all the output. Then the output could be used as a single cell of # a new char.matrix object so one could construct complex tables. x.orig <- x if(!is.matrix(x)) x <- matrix(x, ncol=1) if(!print.it) val <- "" bbox <- my.string.bbox(x) w <- bbox$columns h <- bbox$rows dim(w) <- dim(x) dim(h) <- dim(x) w <- pmax(apply(w, 2, max), min.colwidth) # prefix.width is number of spaces to leave for row names if(missing(prefix.width)) prefix.width <- if(!is.null(dimnames(x)[[1]]) ) max(min.colwidth, min(na.rm = T, median(w), max(nchar( dimnames(x)[[1]])))) else 0 # # If matrix is wider than page split it into 2 pieces, the first of # which is narrow enough, and call print.char.matrix recursively to # print the two pieces. (If second piece is still too wide, further # recursion will take care of it.) # We do not attempt to size to page if print.it=F. # columns.need <- prefix.width + (ncol(x) + 1) * nchar(hsep) + sum(w) if(print.it && (columns.need > page.width)) { which <- (prefix.width + cumsum(w + nchar(hsep)) + nchar(hsep)) <= page.width if(!any(which)) { warning( "page.width too small to print even one column of table" ) which[1] <- T } else if(all(which)) { warning("Internal error in pretty.print") } print(x[, which, drop = F], ..., prefix.width = prefix.width, hsep = hsep, vsep = vsep, csep = csep, print.it = print.it, min.colwidth = min.colwidth, abbreviate.dimnames = abbreviate.dimnames, page.width = columns.need) if(any(!which)) { cat("\n") print(x[, !which, drop = F], ..., prefix.width = prefix.width, hsep = hsep, vsep = vsep, csep = csep, print.it = print.it, min.colwidth = min.colwidth, abbreviate.dimnames = abbreviate.dimnames, page.width = page.width) } } else { # # vsep line is line to print between rows of matrix. It looks like # ------+-----+-----+, made up of vsep and csep characters. # nhsep <- nchar(hsep) if(nchar(vsep) > 0) { if(nchar(csep) != nhsep) stop("nchar(csep) != nhsep") vsep.line <- vector("list", nchar(vsep)) for(j in seq(length = nchar(vsep))) { vsep.line[[j]] <- rep(rep(substring(vsep, j, j), length(w)), w + nhsep) for(i in seq(length = nhsep)) vsep.line[[j]][cumsum(w + nhsep) + i - nhsep] <- substring(csep, i, i) vsep.line[[j]] <- paste(paste(rep(substring( vsep, j, j), prefix.width), collapse = ""), csep, paste(vsep.line[[j]], collapse = ""), "\n", sep = "") } vsep.line <- do.call("paste", c(vsep.line, sep = "", collapse = "")) } else vsep.line <- "" justify <- function(string, width, chop = T) { # # left justify string in a space width characters wide. # string and width may be vectors. # if(chop) string <- substring(string, 1, width) if(nchar(spaces <- " ") < max(width)) spaces <- paste(collapse = "", rep(" ", max( width))) paste(sep = "", string, substring(spaces, 1, width - nchar(string))) } if(!is.null(ndn <- names(dimnames(x)))) { # # print names of dimensions above the matrix # out1 <- paste(sep = "", justify(ndn[1], prefix.width, chop = F), hsep, ndn[2], "\n") if(print.it) cat(out1) else val <- paste(val, out1, sep = "") } # # rn is row names, cn is column names # We must either abbreviate or truncate dimnames to fit in space provided # We might choose to not abbreviate them if we've already abbreviated # them or know that we can fit them into prefix.width spaces. # if(!is.null(rn <- dimnames(x)[[1]])) { if(abbreviate.dimnames && max(nchar(rn)) > prefix.width ) rn <- abbreviate(rn, prefix.width) } else rn <- character(0) if(length(cn <- dimnames(x)[[2]])) { if(abbreviate.dimnames && max(nchar(cn)) > min(w)) cn <- abbreviate(cn, min(w)) out1 <- paste(justify("", prefix.width), hsep, sep = "", paste(collapse = "", sep = "", justify(cn, w), hsep)) out1 <- paste(out1, "\n", sep = "") if(print.it) cat(out1) else val <- paste(val, out1, sep = "") } if(print.it) cat(vsep.line) else val <- paste(val, vsep.line, sep = "") for(i in 1:nrow(x)) { # # loop over rows of matrix # string.break.line breaks strings with newlines into multiple strings # with one line per string. It returns a list with component i containing # the text lines from the i'th input string. # lines <- my.string.break.line(x[i, ]) for(j in 1:max(h[i, ])) { # # loop over text lines in this row # print row name on first line only (rn[i][j] is rn[i] if j==1, "" otherwise) # out1 <- paste(justify(rn[i][j], prefix.width), hsep, sep = "") # linesj is vector of text lines for columns in this line linesj <- unlist(lapply(lines, function(lines, j) lines[j], j)) out1 <- paste(sep = "", out1, paste(sep = "", collapse = hsep, justify(linesj, w, chop = F) ), hsep, "\n") if(print.it) cat(out1) else val <- paste(val, out1, sep = "") } if(print.it) cat(vsep.line) else val <- paste(val, vsep.line, sep = "") } } if(print.it) invisible(x.orig) else val } # # ----------------------------------------------------------------------- # my.string.bbox _ function(s) { # # # Author: S-Plus version 3.2 function # Modified by: Nicky Best # # - returns a 2 element list containing the maximum column width for # each cell in s, and the number of rows in each cell (i.e. 1 + the number # of newline ("\n") symbols in each cell. # # This function mimicks the S-Plus version 3.2 function called # "string.bounding.box", which is used by the "print.char.matrix" function. # The original string.bounding.box function contians a call to a C function # - this has been re-written in S-Plus to ensure that CODA may be run # using earlier versions of S-Plus which do not include this fucntion. # if(!is.matrix(s)) S <- matrix(s, ncol=1) mode(s) <- "character" n <- length(s) out.s <- list() out.s[[1]] <- numeric(ncol(s)) out.s[[2]] <- numeric(nrow(s)) names(out.s) <- c("columns", "rows") k <- 0 for(i in 1:ncol(s)) { for(j in 1:nrow(s)) { k <- k+1 this.s <- substring(s[j, i], first=1:nchar(s[j, i]), last=1:nchar(s[j, i])) out.s$rows[k] <- length(this.s[this.s == "\n"]) + 1 if(out.s$rows[k] > 1) { out.s$columns[k] <- max(diff(c(0, (1:length(this.s))[this.s == "\n"], length(this.s) + 1))) - 1 } else { out.s$columns[k] <- nchar(s[j, i]) } } } out.s } # # ----------------------------------------------------------------------- # my.string.break.line _ function(s) # # Author: S-Plus version 3.2 function # Modified by: Nicky Best # # - returns a list with each element containing a character vector # corresponding to the separate lines of each cell of s # e.g # s <- matrix(c("Hello", "My name\nis Fred"), ncol=2) # my.string.break.line(s) gives: # [[1]]: # [1] "Hello" # # [[2]]: # [1] "My name" "is Fred" # # This function mimicks the S-Plus version 3.2 function called # "string.break.line", which is used by the "print.char.matrix" function. # The original string.break.line function contians a call to a C function # - this has been re-written in S-Plus to ensure that CODA may be run # using earlier versions of S-Plus which do not include this fucntion. # { if(!is.matrix(s)) s <- matrix(s, ncol=1) if(!is.character(s)) s <- as.character(s) n <- length(s) out.s <- list() for(i in 1:n) { this.s <- substring(s[i], first=1:nchar(s[i]), last=1:nchar(s[i])) newline <- c(0, (1:length(this.s))[this.s == "\n"], length(this.s) + 1) for(j in 1:(length(newline) - 1)) { out.s[[i]][j] <- paste(this.s[(newline[j]+1): (newline[j+1]-1)], collapse="") } if(out.s[[i]][length(out.s[[i]])] == "\n") out.s[[i]][length(out.s[[i]])] <- "" } out.s } # # ----------------------------------------------------------------------- # my.format _ function(x, scientific = c(-4, 4), digits = 3, big.mark = "", big.interval = 3, small.mark = "", small.interval = 5, nsmall = 0, justify = "decimal", decimal.mark = ".") # # Author: S-Plus version 3.2 function for formatting vectors for printing # { mod.x <- mode(x) if(mod.x == "missing") x <- NULL if(nargs() == 1 || mod.x != "numeric") { n <- nchar(x) nmax<-max(n) z <- as.character(length(x)) for(i in 1:length(x)) { z[i] <- paste(c(x[i], rep(" ", nmax-n[i])), collapse="") } z } else { expand.small <- function(str, small.mark, small.interval) { if(small.mark == "") { attr(str, "tailsize") <- nchar(str) return(str) } st.nch <- nchar(str) mst <- max(st.nch) # Note next may return str without attr(,"tailsize"). # Is that ok? if(mst <= small.interval) return(str) max.mark <- (mst - 1) %/% small.interval ans <- substring(str, 1, small.interval) for(i in 1:max.mark) { ans <- paste(ans, small.mark, substring(str, i * small.interval + 1, (i + 1) * small.interval), sep = "") } ans <- substring(ans, 1, nchar(ans) - max.mark + ( st.nch - 1) %/% small.interval) attr(ans, "tailsize") <- nchar(ans) ans } charnum <- function(base, expon, decimal.mark, big.mark, big.interval, small.mark, small.interval, nsmall, scientific) { expand.big <- function(str, big.mark, big.interval) { if(big.mark == "") return(str) st.nch <- nchar(str) neg <- substring(str, 1, 1) == "-" mst <- max(st.nch - neg) if(mst <= big.interval) return(str) max.mark <- (mst - 1) %/% big.interval ans <- substring(str, st.nch - big.interval + 1, st.nch) for(i in 1:max.mark) { ans <- paste(substring(str, st.nch - (i + 1) * big.interval + 1, st.nch - (i * big.interval)), big.mark, ans, sep = "") } ans <- substring(ans, 1 + max.mark - (st.nch - 1) %/% big.interval) neg <- substring(ans, 1, 2) == "-," if(any(neg)) ans[neg] <- paste("-", substring(ans[neg], 3), sep = "") ans } if((n <- length(base)) != length(expon)) stop("base and expon must be the same length") special <- base < 0 & expon < 0 base[special] <- abs(base[special]) cbase <- as.character(base) top <- substring(cbase, 1, 2) tail <- substring(cbase, 3) ptop <- substring(top, 2, 2) == "." ptail <- substring(tail, 1, 1) == "." ntop <- substring(top, 1, 1 + (!ptop)) ntail <- substring(tail, 1 + ptail) ans <- final.top <- final.tail <- character(n) zero <- expon == 0 if(any(zero)) { final.tail[zero] <- ntail[zero] final.top[zero] <- ntop[zero] } pos <- expon > 0 & expon < nchar(ntail) if(any(pos)) { final.top[pos] <- paste(ntop[pos], substring( ntail[pos], 1, expon[pos]), sep = "") final.tail[pos] <- substring(ntail[pos], expon[ pos] + 1) } for(i in (1:n)[!(pos | zero)]) { if(expon[i] < 0) { final.top[i] <- "0" final.tail[i] <- paste(paste(rep("0", abs( expon[i]) - 1), collapse = ""), ntop[i], ntail[i], sep = "") } else { final.top[i] <- paste(ntop[i], ntail[i], paste(rep("0", expon[i] - nchar(ntail[i])), collapse = ""), sep = "") final.tail[i] <- "" } } if(!is.null(nsmall) && nsmall >= 1) { nsmall <- round(nsmall) tail.nc <- nchar(final.tail) if(any(tail.nc < nsmall)) final.tail <- substring(paste(final.tail, rep( paste(rep("0", nsmall), collapse = ""), n), sep = ""), 1, nsmall) } tmp <- get("expand.small", frame = sys.parent())( final.tail, small.interval = small.interval, small.mark = small.mark) ntmp <- nchar(tmp) if(any(ntmp == 0)) ans[ntmp == 0] <- expand.big(final.top[ntmp == 0], big.mark = big.mark, big.interval = big.interval) if(any(ntmp != 0)) ans[ntmp != 0] <- paste(expand.big(final.top[ ntmp != 0], big.mark = big.mark, big.interval = big.interval), decimal.mark, tmp[ntmp != 0 ], sep = "") ans[special] <- paste("-", ans[special], sep = "") attr(ans, "tailsize") <- ntmp ans } # # start of main function # if(non.fin <- any(non.fin.val <- !is.finite(x))) { big.ans <- as.character(x) x <- x[!non.fin.val] if(is.logical(scientific) && length(scientific) > 1) { scientific <- rep(scientific, length = length( big.ans)) scientific <- scientific[!non.fin.val] } } n <- length(x) if(!n) { if(non.fin) return(big.ans) else return(character(0)) } adj.ind <- charmatch(justify, c("none", "left", "right", "decimal"), nomatch = NA) if(is.na(adj.ind)) stop(paste("unknown value for justify:", justify)) else justify <- c("none", "left", "right", "decimal")[adj.ind] x.att <- attributes(x) x <- as.vector(x) tailsize <- numeric(n) expon <- floor(log10(abs(x))) expon[x == 0] <- 0 basenum <- round(x * 10^( - expon), digits - 1) if(any(is.infinite(basenum))) { si <- sum(is.infinite(basenum)) if(si == 1) stop("one number out of range") else stop(paste(si, "numbers out of range")) } # now do adjustment because rounding may have changed something expon <- ifelse(abs(basenum) >= 10, expon + 1, ifelse(abs( basenum) < 1 & basenum != 0, expon - 1, expon)) basenum <- ifelse(abs(basenum) >= 10, basenum/10, ifelse(abs( basenum) < 1, basenum * 10, basenum)) if(is.logical(scientific)) science.use <- rep(scientific, length = n) else science.use <- expon >= scientific[2] | expon <= scientific[1] ans <- character(n) if(any(science.use)) { neg <- basenum[science.use] < 0 if(!is.null(nsmall) && nsmall >= 1) { tmp <- substring(basenum[science.use], 3 + neg) nsmall <- round(nsmall) tail.nc <- nchar(tmp) if(any(tail.nc < nsmall)) { tmp <- substring(paste(tmp, rep(paste(rep("0", nsmall), collapse = ""), sum(science.use)), sep = ""), 1, nsmall) attr(tmp, "tailsize") <- nchar(tmp) } else { tmp <- expand.small(tmp, small.interval = small.interval, small.mark = small.mark) attr(tmp, "tailsize") <- nchar(tmp) } } else { tmp <- expand.small(substring(basenum[ science.use], 3 + neg), small.interval = small.interval, small.mark = small.mark) } tailsize[science.use] <- attr(tmp, "tailsize") + 1 + nchar(paste(expon[science.use])) ans[science.use] <- paste(substring(basenum[science.use ], 1, neg + 1), decimal.mark, tmp, "e", expon[ science.use], sep = "") } if(any(!science.use)) { tmp <- charnum(basenum[!science.use], expon[! science.use], decimal.mark = decimal.mark, big.interval = big.interval, big.mark = big.mark, small.interval = small.interval, small.mark = small.mark, nsmall = nsmall, scientific = scientific) tailsize[!science.use] <- attr(tmp, "tailsize") ans[!science.use] <- tmp } ans <- switch(justify, none = ans, right = { tailsize <- nchar(ans) ndiff <- diff(range(tailsize)) substring(paste(rep(paste(rep(" ", ndiff), collapse = ""), n), ans, sep = ""), tailsize - min(tailsize) + 1) } , left = { tailsize <- nchar(ans) ndiff <- diff(range(tailsize)) substring(paste(ans, rep(paste(rep(" ", ndiff), collapse = ""), n), sep = ""), 1, max( tailsize)) } , decimal = { nc.str <- nchar(ans) ndiff <- diff(range(tailsize)) if(ndiff) ans <- substring(paste(ans, rep(paste(rep(" ", ndiff + 1), collapse = ""), n), sep = ""), 1, nc.str + max(tailsize) - tailsize + ( tailsize == 0)) topsize <- nc.str - tailsize - (tailsize != 0) ndiff <- diff(range(topsize)) if(ndiff) ans <- substring(paste(rep(paste(rep(" ", ndiff), collapse = ""), n), ans, sep = ""), topsize - min(topsize) + 1) ans } ) if(non.fin) { big.ans[!non.fin.val] <- ans attributes(big.ans) <- x.att big.ans } else { attributes(ans) <- x.att ans } } } # # # ----------------------------------------------------------------------- # working.data _ function() { # # working.data -- Set up working data array defined by current data.defaults(): # Array has 3 dimensions corresponding to: # [1] iterations # [2] chains # [3] parameters # Also, create object `parmsmat.iter' giving iteration # numbers currently selected, and `parmnames' - an object # of class "char.matrix" containing variable labels to be # used when printing any tables # # Author: Nicky Best # # # It would be possible to use the <<- method of assigning values # to permanent objects stored in the working database. # However, all such assignments remain pending until the completion # of the top level expression (which in this case is the CODA() # function itself). Hence S-Plus will retain a copy of each previous # value of objects such as parmsmat - this causes excessive memory # buildup inside S-Plus, eventually causing CODA to crash if the .out # files are large. # This mechanism can be overridden by using `assign' with the # `immediate=T' option (see p. 4-25 of the S-Plus version 3.2 # Programmers Manual). # if(any(c(changes()$parms, changes()$chains, changes()$iter, changes()$thin))) { if(!globals()$first.time) { cat("\nRe-creating working data object....\n") } iters <- seq(data.defaults()$start, data.defaults()$end, data.defaults()$thin) globals(Niters = length(iters)) globals(Nchains = length(data.defaults()$chains)) globals(Nparms = length(data.defaults()$parms)) tmp.expr <- "array(c(" for(i in 1:globals()$Nchains) { tmp.expr <- paste(tmp.expr, paste("globals()$out[[", data.defaults()$chains[i], "]][iters, data.defaults()$parms+1]", sep=""), sep=",") } tmp.expr <- paste(tmp.expr, "), dim=c(", globals()$Niters, ",", globals()$Nparms, ", ", globals()$Nchains, "))", sep="") globals(parmsmat = eval(parse(text=tmp.expr))) # # Create new vector giving iteration numbers # corresponding to rows of parmsmat # globals(parmsmat.iter = iters+globals()$out[[1]][1, "iter"]) globals(parmnames = pretty.print(matrix(c("VARIABLE", "========", "", globals()$Labels[data.defaults()$parms]), ncol=1), print.it=F, hsep="", vsep="",csep=" ")) globals(parmnames2 = pretty.print(matrix(c("", "VARIABLE", "========", "", globals()$Labels[data.defaults()$parms]), ncol=1), print.it=F, hsep="", vsep="",csep=" ")) } # # Create new indicator for positive parameters in parmsmat # NB: this needs to be done if either parmsmat has been changed # or parmsmat remains the same but the positive variables have changed # if(length(func.defaults()$pos.parms)==0) { func.defaults(pos.parmsmat = 0) } else { if(func.defaults()$pos.parms[1]==0) { func.defaults(pos.parmsmat = 0) } else { func.defaults(pos.parmsmat = numeric()) for(i in 1:globals()$Nparms) { if(any(func.defaults()$pos.parms==data.defaults()$parms[i])) func.defaults(pos.parmsmat = append(func.defaults()$pos.parmsmat, i)) } } } if(length(func.defaults()$pos.parmsmat)==0) func.defaults(pos.parmsmat = 0) # # Create new indicator for (0, 1) parameters in parmsmat # if(length(func.defaults()$prop.parms)==0) { func.defaults(prop.parmsmat = 0) } else { if(func.defaults()$prop.parms[1]==0) { func.defaults(prop.parmsmat = 0) } else { func.defaults(prop.parmsmat = numeric()) for(i in 1:globals()$Nparms) { if(any(func.defaults()$prop.parms == data.defaults()$parms[i])) func.defaults(prop.parmsmat = append(func.defaults()$prop.parmsmat, i)) } } } if(length(func.defaults()$prop.parmsmat)==0) func.defaults(prop.parmsmat = 0) changes(parms = F, chains = F, iter = F, thin = F, pos.parms = F) } # # ----------------------------------------------------------------------- # open.graphics _ function() { # # open.graphics -- Opens a graphics window selected by the user according # to the operating system being used to run S-Plus # # Author: Nicky Best # restart(T) if(!globals()$graph.open) { cat("\nSelect graphics device required:\n") graphics.menu <- c("openlook (UNIX)", "motif (UNIX)", "X11 (UNIX)", "win.graph (DOS)") globals(which.dev = menu(graphics.menu)) cat("\nOpening graphics window...\n") switch(globals()$which.dev, # # UNIX options: # ============ # openlook(), motif(), X11(), # # PC option: # ========== # win.graph() ) globals(graph.open = T) } } # # ----------------------------------------------------------------------- # plot.to.file _ function() { # # plot.to.file -- gives user the option to save plots currently appearing # in graphics window to a postscipt file # # Author: Nicky Best # # Modified by NGB 8/4/97 to enable graphics to be saved to file when # using S-Plus for Windows. Thanks to Siem Heisterkamp for suggesting # this fix. # if(globals()$which.dev < 4) { while(T) { # UNIX cat("\nDo you want to save current plots as a postscript file (y/n) ?\n") ans <- scan(what = character(), n = 1, strip.white=T) if(length(ans) > 0) break } } else { # WINDOWS cat("\nDo you want to save current plots to a file?\n") plot.to.file.menu <- c("No", "Yes - as a Placeable Windows Metafile", "Yes - as a Postscript file") which.file <- menu(plot.to.file.menu) if(which.file == 1) ans <- "n" else ans <- "y" } if(ans == "Y" | ans == "y") { while(T) { cat("\nEnter name you want to call this file:\n") file.name <- scan(what = character(), n = 1, strip.white=T) if(length(file.name) > 0) { file.name <- paste(globals()$Home, file.name, sep="") break } } if(globals()$which.dev < 4) { # # UNIX graphics device # tmp <- switch(func.defaults()$ps.plot, printgraph(method = "postscript", file = file.name, horizontal=F, print=F, paper="a4", width=6, height=9,maximize=F), printgraph(method = "postscript", file = file.name, horizontal=T, print=F, paper="a4", width=9, height=6,maximize=F), printgraph(method = "postscript", file = file.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) } } else { # # WINDOWS graphics device # tmp1 <- switch(which.file-1, { # metafile dev.copy(device = win.printer, file = file.name, format = "placeable metafile", width=9, height=6) dev.off() }, { # postscript file tmp3<-switch(func.defaults()$ps.plot, dev.copy(postscript, file = file.name, horizontal=F, width=6, height=9), dev.copy(postscript, file = file.name, horizontal=T, width=9, height=6), dev.copy(postscript, file = file.name, horizontal=F, width=6, height=4) ) if(is.null(tmp3)) { # default option dev.copy(postscript, file = file.name, horizontal=F, width=6, height=9) } dev.off() } ) } } } # # ----------------------------------------------------------------------- # pause _ function(object, ...) { # # pause -- function to create a pause between multiple pages of plots # to enable user to view/print plots before next page appears # # Author: Nicky Best # # This is basically the S-Plus `browser()' function, but with # different prompt and message options. An option to print # graphics object to a postscript file is also included # tmp <- plot.to.file() if(nargs()) UseMethod("browser") else { nframe <- sys.parent() msg <- paste(deparse(sys.call(nframe)), collapse = " ") browser.default(nframe, message = "\nPlots span more than 1 page", prompt = "Type 0 to continue scrolling: ") } } # # ----------------------------------------------------------------------- # 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(!is.null(globals()$out)) { while(T) { cat("\nDo you want to save the BUGS output as an S-Plus object file(y/n) ?\n") ans <- scan(what = character(), n = 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(), n = 1, strip.white=T) assign(fname, globals()$out, where=1) } } graphics.off() globals.remove() cat("Have a BUGS-free day!\n") }