Rev 5661 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
cumuplot <- function(x, probs=c(0.025,0.5,0.975), ylab="", lty=c(2,1),lwd=c(1,2), type="l", ask,auto.layout=TRUE, col=1, ...){if (missing(ask)) {ask <- if (is.R()) {dev.interactive()}else {interactive()}}cquantile <- function(z, probs){## Calculates cumulative quantile of a vectorcquant <- matrix(0, nrow=length(z), length(probs))for(i in seq(along=z)) # for loop proved faster than apply hereif (is.R()) {cquant[i,] <- quantile(z[1:i], probs=probs, names=FALSE)}else{cquant[i,] <- quantile(z[1:i], probs=probs)}cquant <- as.data.frame(cquant)names(cquant) <- paste(formatC(100*probs,format="fg",width=1,digits=7),"%", sep="") # just like quantile.defaultreturn(cquant)}oldpar <- NULLon.exit(par(oldpar))if (auto.layout) {oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x),Nparms = nvar(x)))}if (!is.mcmc.list(x))x <- mcmc.list(as.mcmc(x))Iterations <- as.vector(time(x))for (i in 1:nchain(x)) {for (j in 1:nvar(x)) {Y <- cquantile(as.matrix(x[[i]])[,j], probs=probs)if (!is.R()) Y <- as.matrix(Y)matplot(Iterations, Y, ylab=ylab, lty=lty, lwd=lwd, type=type,col=col, ...)title(paste(varnames(x)[j], ifelse(is.null(chanames(x)),"", ":"), chanames(x)[i], sep = ""))if (i == 1 & j == 1)oldpar <- c(oldpar, par(ask=ask))}}}