1"[.mcmc" <- function (x, i, j, drop = missing(i)) 2{ 3 ## In S-PLUS the code is altered so that the user can 4 ## pick out particular parameters by calling mcmc.obj[,c("param1", "param2")] 5 xstart <- start(x) 6 xthin <- thin(x) 7 if (is.R()) { 8 y <- NextMethod("[") 9 } 10 else { 11 y <- as.matrix(x)[i,j] 12 } 13 if (length(y) == 0 || is.null(y)) 14 return(y) 15 if (missing(i)) 16 return(mcmc(y, start = xstart, thin = xthin)) 17 else 18 return(y) 19} 20 21"as.mcmc" <- function (x, ...) 22 UseMethod("as.mcmc") 23 24"as.mcmc.default" <- function (x, ...) 25 if (is.mcmc(x)) x else mcmc(x) 26 27"as.ts.mcmc" <- function (x, ...) 28{ 29 x <- as.mcmc(x) 30 y <- ts(x, start = start(x), end = end(x), deltat = thin(x)) 31 attr(y, "mcpar") <- NULL 32 return(y) 33} 34 35"start.mcmc" <- function (x, ...) 36{ 37 mcpar(as.mcmc(x))[1] 38} 39 40"end.mcmc" <- function (x, ...) 41{ 42 mcpar(as.mcmc(x))[2] 43} 44 45"frequency.mcmc" <- function (x, ...) 46{ 47 1/thin.mcmc(x) 48} 49 50"thin.mcmc" <- function (x, ...) 51{ 52 mcpar(as.mcmc(x))[3] 53} 54 55"is.mcmc" <- function (x) 56{ 57 if (inherits(x, "mcmc")) 58 if (length(dim(x)) == 3) 59 stop("Obsolete mcmc object\nUpdate with a command like\nx <- mcmcUpgrade(x)") 60 else TRUE 61 else FALSE 62} 63 64"mcmc" <- function (data = NA, start = 1, end = numeric(0), thin = 1) 65{ 66 if (is.matrix(data)) { 67 niter <- nrow(data) 68 nvar <- ncol(data) 69 } 70 else if (is.data.frame(data)) { 71 if (!all(sapply(data, is.numeric))) { 72 stop ("Data frame contains non-numeric values") 73 } 74 data <- as.matrix(data) 75 niter <- nrow(data) 76 nvar <- ncol(data) 77 } 78 else { 79 niter <- length(data) 80 nvar <- 1 81 } 82 thin <- round(thin) 83 if (length(start) > 1) 84 stop("Invalid start") 85 if (length(end) > 1) 86 stop("Invalid end") 87 if (length(thin) != 1) 88 stop("Invalid thin") 89 if (missing(end)) 90 end <- start + (niter - 1) * thin 91 else if (missing(start)) 92 start <- end - (niter - 1) * thin 93 nobs <- floor((end - start)/thin + 1.0) ### patch 94 if (niter < nobs) 95 stop("Start, end and thin incompatible with data") 96 else { 97 end <- start + thin * (nobs - 1) 98 if (nobs < niter) 99 data <- data[1:nobs, , drop = FALSE] 100 } 101 attr(data, "mcpar") <- c(start, end, thin) 102 attr(data, "class") <- "mcmc" 103 data 104} 105 106"print.mcmc" <- function (x, ...) 107{ 108 x.orig <- x 109 cat("Markov Chain Monte Carlo (MCMC) output:\nStart =", start(x), 110 "\nEnd =", end(x), "\nThinning interval =", thin(x), "\n") 111 attr(x, "mcpar") <- NULL 112 attr(x, "class") <- NULL 113 NextMethod("print", ...) 114 invisible(x.orig) 115} 116 117 118"as.matrix.mcmc" <- function (x, iters = FALSE, ...) 119{ 120 y <- matrix(nrow = niter(x), ncol = nvar(x) + iters) 121 var.cols <- iters + 1:nvar(x) 122 if (iters) 123 y[, 1] <- as.vector(time(x)) 124 y[, var.cols] <- x 125 rownames <- character(ncol(y)) 126 if (iters) 127 rownames[1] <- "ITER" 128 rownames[var.cols] <- varnames(x, allow.null = FALSE) 129 dimnames(y) <- list(NULL, rownames) 130 return(y) 131} 132 133"time.mcmc" <- function (x, ...) 134{ 135 x <- as.mcmc(x) 136 ts(seq(from = start(x), to = end(x), by = thin(x)), start = start(x), 137 end = end(x), deltat = thin(x)) 138} 139 140"window.mcmc" <- function (x, start, end, thin, ...) 141{ 142 ts.eps <- getOption("ts.eps") 143 xmcpar <- mcpar(x) 144 xstart <- xmcpar[1] 145 xend <- xmcpar[2] 146 xthin <- xmcpar[3] 147 if (missing(thin)) 148 thin <- xthin 149 else if (thin%%xthin != 0) { 150 thin <- xthin 151 warning("Thin value not changed") 152 } 153 xtime <- as.vector(time(x)) 154 if (missing(start)) 155 start <- xstart 156 else if (length(start) != 1) 157 stop("bad value for start") 158 else if (start < xstart) { 159 start <- xstart 160 warning("start value not changed") 161 } 162 if (missing(end)) 163 end <- xend 164 else if (length(end) != 1) 165 stop("bad value for end") 166 else if (end > xend) { 167 end <- xend 168 warning("end value not changed") 169 } 170 if (start > end) 171 stop("start cannot be after end") 172 if (all(abs(xtime - start) > abs(start) * ts.eps)) { 173 start <- xtime[(xtime > start) & ((start + xthin) > xtime)] 174 } 175 if (all(abs(end - xtime) > abs(end) * ts.eps)) { 176 end <- xtime[(xtime < end) & ((end - xthin) < xtime)] 177 } 178 use <- 1:niter(x) 179 use <- use[use >= trunc((start - xstart)/xthin + 1.5) & 180 use <= trunc((end - xstart)/xthin + 1.5) & 181 (use - trunc((start- xstart)/xthin + 1.5))%%(thin%/%xthin) == 0] 182 y <- if (is.matrix(x)) 183 x[use, , drop = FALSE] 184 else x[use] 185 return(mcmc(y, start=start, end=end, thin=thin)) 186} 187 188"mcpar" <- function (x) 189{ 190 attr(x, "mcpar") 191} 192 193"mcmcUpgrade" <- function (x) 194{ 195 if (inherits(x, "mcmc")) { 196 if (length(dim(x)) == 3) { 197 nchain <- dim(x)[3] 198 xtspar <- attr(x, "tspar") 199 xstart <- xtspar[1] 200 xend <- xtspar[2] 201 xthin <- xtspar[3] 202 out <- vector("list", nchain) 203 for (i in 1:nchain) { 204 y <- unclass(x)[, , 1, drop = TRUE] 205 attr(y, "title") <- NULL 206 attr(y, "tspar") <- NULL 207 out[[i]] <- mcmc(y, start = xstart, end = xend, 208 thin = xthin) 209 } 210 if (nchain == 1) 211 return(out[[1]]) 212 else return(mcmc.list(out)) 213 } 214 else return(x) 215 } 216 else stop("Can't upgrade") 217} 218 219"thin" <- 220function (x, ...) 221 UseMethod("thin") 222 223"set.mfrow" <- 224function (Nchains = 1, Nparms = 1, nplots = 1, sepplot = FALSE) 225{ 226 ## Set up dimensions of graphics window: 227 ## If only density plots OR trace plots are requested, dimensions are: 228 ## 1 x 1 if Nparms = 1 229 ## 1 X 2 if Nparms = 2 230 ## 2 X 2 if Nparms = 3 or 4 231 ## 3 X 2 if Nparms = 5 or 6 or 10 - 12 232 ## 3 X 3 if Nparms = 7 - 9 or >= 13 233 ## If both density plots AND trace plots are requested, dimensions are: 234 ## 1 x 2 if Nparms = 1 235 ## 2 X 2 if Nparms = 2 236 ## 3 X 2 if Nparms = 3, 5, 6, 10, 11, or 12 237 ## 4 x 2 if Nparms otherwise 238 ## If separate plots are requested for each chain, dimensions are: 239 ## 1 x 2 if Nparms = 1 & Nchains = 2 240 ## 2 X 2 if Nparms = 2 & Nchains = 2 OR Nparms = 1 & Nchains = 3 or 4 241 ## 3 x 2 if Nparms = 3 or >= 5 & Nchains = 2 242 ## OR Nchains = 5 or 6 or 10 - 12 (and any Nparms) 243 ## 2 x 3 if Nparms = 2 or 4 & Nchains = 3 244 ## 4 x 2 if Nparms = 4 & Nchains = 2 245 ## OR Nchains = 4 & Nparms > 1 246 ## 3 x 3 if Nparms = 3 or >= 5 & Nchains = 3 247 ## OR Nchains = 7 - 9 or >= 13 (and any Nparms) 248 mfrow <- if (sepplot && Nchains > 1 && nplots == 1) { 249 ## Separate plots per chain 250 ## Only one plot per variable 251 if (Nchains == 2) { 252 switch(min(Nparms, 5), 253 c(1,2), 254 c(2,2), 255 c(3,2), 256 c(4,2), 257 c(3,2)) 258 } 259 else if (Nchains == 3) { 260 switch(min(Nparms, 5), 261 c(2,2), 262 c(2,3), 263 c(3,3), 264 c(2,3), 265 c(3,3)) 266 } 267 else if (Nchains == 4) { 268 if (Nparms == 1) 269 c(2,2) 270 else 271 c(4,2) 272 } 273 else if (any(Nchains == c(5,6,10,11,12))) 274 c(3,2) 275 else if (any(Nchains == c(7,8,9)) || Nchains >=13) 276 c(3,3) 277 278 } 279 else { 280 if (nplots==1) { 281 ## One plot per variable 282 mfrow <- switch(min(Nparms,13), 283 c(1,1), 284 c(1,2), 285 c(2,2), 286 c(2,2), 287 c(3,2), 288 c(3,2), 289 c(3,3), 290 c(3,3), 291 c(3,3), 292 c(3,2), 293 c(3,2), 294 c(3,2), 295 c(3,3)) 296 } 297 else { 298 ## Two plot per variable 299 ## 300 mfrow <- switch(min(Nparms, 13), 301 c(1,2), 302 c(2,2), 303 c(3,2), 304 c(4,2), 305 c(3,2), 306 c(3,2), 307 c(4,2), 308 c(4,2), 309 c(4,2), 310 c(3,2), 311 c(3,2), 312 c(3,2), 313 c(4,2)) 314 } 315 } 316 return(mfrow) 317} 318 319head.mcmc <- function(x, n = 6L, ...) { 320 window.mcmc(x, end=min(start.mcmc(x) + n * thin.mcmc(x), end.mcmc(x))) 321} 322 323tail.mcmc <- function(x, n = 6L, ...) { 324 window.mcmc(x, start=max(end.mcmc(x) - n * thin.mcmc(x), start.mcmc(x))) 325} 326