### This file consolidates all of the summary functions supplied by ASP. ## Summaries ##This is the usual collection of coda functions I like to apply to an ##model. "grandSummary" <- function(param) { if (all(is.na(param))) { cat("Data not available for this parameter set.\n") } else { print(try(summary(param))) print(try(gelman.diag(param,transform=TRUE))) cat ("Autocorrelations\n") print(try(autocorr.diag(param))) cat ("Effective Sizes\n") print(try(effectiveSize(param))) cat ("RejectionRate\n") print(try(rejectionRate(param))) ## Return a dividing line which will get printed. cat("-------------------------------------------\n") NULL } } safeGelman.diag <- function(param,confidence=.95, transform=TRUE,autoburnin=FALSE) { result <- try(gelman.diag(param,confidence,transform,autoburnin)) if (class(result) == "try-error") result <- NA if (class(result) == "try") result <- NA result } mpsrf <- function(param,confidence=.95,transform=TRUE,autoburnin=FALSE) { result <- safeGelman.diag(param,confidence,transform,autoburnin)$mpsrf if (is.null(result)) { NA } else { Re(result) ## Often comes up as complex number. } } maxpsrf <- function(param,confidence=.95,transform=TRUE,autoburnin=FALSE) { result <- safeGelman.diag(param,confidence,transform,autoburnin)$psrf if (is.null(result)) { NA } else { max(result[,1]) } } safeAutocorr.diag <- function(param,lags=c(0,1,5,10,50),relative=TRUE) { result <- try(autocorr.diag.mcmc.list(param,lags=lags,relative=relative)) if (class(result) == "try-error") result <- NA if (class(result) == "try") result <- NA result } safeEffectiveSize <- function (x) { result <- try(effectiveSize(x)) if (class(result) == "try-error") result <- NA if (class(result) == "try") result <- NA result } safeRejectionRate <- function (x) { result <- try(rejectionRate(x)) if (class(result) == "try-error") result <- NA if (class(result) == "try") result <- NA result } meanVar <- function (obj,pnames=rownames(stats)) { if (all(is.na(obj))) { NA } else { stats <- summary(obj)$statistics[,1:2] snames <- colnames(stats) result <- as.vector(t(stats)) names(result) <- paste(rep(pnames,each=2),rep(snames,length(snames)),sep=".") result } } mvr <- function (obj,pnames) { if (length(obj)==1 || class(obj)[1]=="asp.table") { rep(NA,2*length(pnames)+1) } else { if (nvar(obj) < length(pnames)) { mv <- meanVar(obj,pnames[1:nvar(obj)]) mv <- c(mv,rep(NA,2*(length(pnames)-nvar(obj)))) } else if (nvar(obj) > length(pnames)) { mv <- meanVar(obj[,1:length(pnames)],pnames) } else { mv <- meanVar(obj,pnames) } R <- safeGelman.diag(obj)$mpsrf if (is.null(R)) { R <- NA } else { R <- Re(R) # For some reason this is often complex. } c(mv,R=R) } } RNC <- function (ps) { R <- mpsrf(ps,transform=TRUE,autoburnin=FALSE) N <- min(effectiveSize(ps)) C <- max(autocorr.diag(ps,lags=5)) list(R=R,N=N,C=C) } maxRNC <- function(pslist) { rnclist <- lapply(pslist,RNC) R <- do.call("max",lapply(rnclist,function(x) x$R)) N <- do.call("min",lapply(rnclist,function(x) x$N)) C <- do.call("max",lapply(rnclist,function(x) x$C)) list(R=R,N=N,C=C) }