File: cumuplot.R

package info (click to toggle)
r-cran-coda 0.13-2-1
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 456 kB
  • sloc: makefile: 2
file content (54 lines) | stat: -rw-r--r-- 1,698 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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 vector
        cquant <- matrix(0, nrow=length(z), length(probs))
        for(i in seq(along=z))  # for loop proved faster than apply here
            if (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",wid=1,digits=7),
                               "%", sep="")  # just like quantile.default
        return(cquant)
    }

    oldpar <- NULL
    on.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 <- 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))
        }
    }
}