1 2 3### Copyright (C) 2005 Deepayan Sarkar 4### <Deepayan.Sarkar@R-project.org>, Douglas Bates 5### <Douglas.Bates@R-project.org>. See file COPYING for license 6### terms. 7 8 9 10 11 12### unexported helper function to obtain a valid subset argument. 13### Mostly to ensure that we don't up with a huge vector of integer 14### indices in the trivial cases. 15 16 17 18thinned.indices <- function(object, n = NROW(object), start = 1, thin = 1) 19{ 20 if (is.mcmc(object) && 21 (start * thin != 1) && 22 !all(mcpar(object)[-2] == 1)) 23 warning("mcmc object is already thinned") 24 if (start < 1) stop("Invalid start") 25 else if (start == 1) 26 { 27 if (thin < 1) stop("Invalid thin") 28 else if (thin == 1) TRUE 29 else rep(c(TRUE, FALSE), c(1, thin-1)) 30 } 31 else 32 { 33 if (thin < 1) stop("Invalid thin") 34 else if (thin == 1) -seq(length = start-1) 35 else start + thin * (0:(floor(n - start) / thin)) 36 } 37} 38 39 40 41## most functions will have methods for mcmc and mcmclist objects. In 42## the second case, another grouping variable is added. By default, 43## this will be used for ``grouped displays'' (when possible), but 44## could also be used for conditioning by setting groups=FALSE 45 46 47 48 49 50## levelplot, analog of plot.crosscorr. Defaults changed to match 51## those of plot.crosscorr as much as possible. 52 53 54levelplot.mcmc <- 55 function(x, data = NULL, 56 main = attr(x, "title"), 57 start = 1, thin = 1, 58 ..., 59 xlab = "", ylab = "", 60 cuts = 10, 61 at = do.breaks(c(-1.001, 1.001), cuts), 62 col.regions = topo.colors(100), 63 subset = thinned.indices(x, start = start, thin = thin)) 64{ 65 if (!is.R()) { 66 stop("This function is not yet available in S-PLUS") 67 } 68 cormat <- cor(x[subset, ]) 69 cormat <- cormat[, rev(seq(length = ncol(cormat)))] 70 levelplot(cormat, 71 main = main, ..., 72 cuts = cuts, at = at, 73 col.regions = col.regions, 74 xlab = xlab, ylab = ylab) 75} 76 77 78## the mcmc.list method wouldn't do any grouping (so should have 79## outer=TRUE by default). It hasn't been written yet because 80 81 82 83 84## The splom (FIXME: not yet written) method may be more useful 85 86## in progress, unexported. Planning to make it be like levelplot on 87## the lower diagonal (maybe ellipses instead of plain boxes) and 88## normal splom on the upper diagonal. Not much point in having tick 89## marks. Names in the middle save space, unlike in levelplot which 90## stupidly shows correlation=1 on the diagonal. 91 92 93splom.mcmc <- 94 function(x, data = NULL, 95 main = attr(x, "title"), 96 start = 1, thin = 1, 97 as.matrix = TRUE, 98 xlab = "", ylab = "", 99 cuts = 10, 100 at = do.breaks(c(-1.001, 1.001), cuts), 101 col.regions = topo.colors(100), 102 ..., 103 pscales = 0, 104 subset = thinned.indices(x, start = start, thin = thin)) 105{ 106## cormat <- cor(x[subset, ]) 107## cormat <- cormat[, rev(seq(length = ncol(cormat)))] 108 if (!is.R()) { 109 stop("This function is not yet available in S-PLUS") 110 } 111 splom(as.data.frame(x[subset, ]), 112 as.matrix = as.matrix, 113 main = main, ..., 114 pscales = pscales, 115 cuts = cuts, at = at, 116 lower.panel = function(x, y, ...) { 117 corval <- cor(x, y) 118 grid::grid.text(lab = round(corval, 2)) 119 }, 120 col.regions = col.regions, 121 xlab = xlab, ylab = ylab) 122} 123 124 125 126 127 128 129 130 131### methods for densityplot (mcmc and mcmc.list) 132 133 134 135densityplot.mcmc <- 136 function(x, data = NULL, 137 outer, aspect = "xy", 138 default.scales = list(relation = "free"), 139 start = 1, thin = 1, 140 main = attr(x, "title"), 141 xlab = "", 142 plot.points = "rug", 143 ..., 144 subset = thinned.indices(x, start = start, thin = thin)) 145{ 146 if (!is.R()) { 147 stop("This function is not yet available in S-PLUS") 148 } 149 if (!missing(outer)) warning("specification of outer ignored") 150 data <- as.data.frame(x) 151 form <- 152 as.formula(paste("~", 153 paste(lapply(names(data), as.name), 154 collapse = "+"))) 155 156 157### The following is one possible approach, but it does not generalize 158### to mcmclist objects. 159 160 161## densityplot(form, data = data, 162## outer = outer, 163## aspect = aspect, 164## default.scales = default.scales, 165## main = main, 166## xlab = xlab, 167## plot.points = plot.points, 168## subset = eval(subset), 169## ...) 170 171 172### This one does, with the only downside I can think of being that 173### subscripts, if used, will give indices in subsetted data, not 174### original. But that's true even if the original mcmc object was 175### itself already thinned. 176 177 densityplot(form, data = data[subset, , drop=FALSE], 178 outer = TRUE, 179 aspect = aspect, 180 default.scales = default.scales, 181 main = main, 182 xlab = xlab, 183 plot.points = plot.points, 184 ...) 185} 186 187 188densityplot.mcmc.list <- 189 function(x, data = NULL, 190 outer = FALSE, groups = !outer, 191 aspect = "xy", 192 default.scales = list(relation = "free"), 193 start = 1, thin = 1, 194 main = attr(x, "title"), 195 xlab = "", 196 plot.points = "rug", 197 ..., 198 subset = thinned.indices(x[[1]], start = start, thin = thin)) 199{ 200 if (!is.R()) { 201 stop("This function is not yet available in S-PLUS") 202 } 203 if (groups && outer) warning("'groups=TRUE' ignored when 'outer=TRUE'") 204 datalist <- lapply(x, function(x) as.data.frame(x)[subset, ,drop=FALSE]) 205 data <- do.call("rbind", datalist) 206 form <- 207 if (outer) 208 as.formula(paste("~", 209 paste(lapply(names(data), as.name), 210 collapse = "+"), 211 "| .run")) 212## as.formula(paste("~", 213## paste(names(data), 214## collapse = "+"), 215## "| .run")) 216 else 217 as.formula(paste("~", 218 paste(lapply(names(data), as.name), 219 collapse = "+"))) 220## as.formula(paste("~", 221## paste(names(data), 222## collapse = "+"))) 223 ##data[["index"]] <- seq(length = nrow(x[[1]]))[subset] 224 .run <- gl(length(datalist), nrow(datalist[[1]])) 225 if (groups && !outer) 226 densityplot(form, data = data, 227 outer = TRUE, 228 groups = .run, 229 aspect = aspect, 230 default.scales = default.scales, 231 main = main, 232 xlab = xlab, 233 plot.points = plot.points, 234 ...) 235 else 236 densityplot(form, data = data, 237 outer = TRUE, 238 aspect = aspect, 239 default.scales = default.scales, 240 main = main, 241 xlab = xlab, 242 plot.points = plot.points, 243 ...) 244} 245 246 247### methods for qqmath (mcmc and mcmc.list) 248 249 250 251 252qqmath.mcmc <- 253 function(x, data = NULL, 254 outer, aspect = "xy", 255 default.scales = list(y = list(relation = "free")), 256 prepanel = prepanel.qqmathline, 257 start = 1, thin = 1, 258 main = attr(x, "title"), 259 ylab = "", 260 ..., 261 subset = thinned.indices(x, start = start, thin = thin)) 262{ 263 if (!is.R()) { 264 stop("This function is not yet available in S-PLUS") 265 } 266 if (!missing(outer)) warning("specification of outer ignored") 267 data <- as.data.frame(x) 268 form <- 269 as.formula(paste("~", 270 paste(lapply(names(data), as.name), 271 collapse = "+"))) 272 qqmath(form, data = data[subset, ,drop=FALSE], 273 outer = TRUE, 274 aspect = aspect, 275 prepanel = prepanel, 276 default.scales = default.scales, 277 main = main, 278 ylab = ylab, 279 ...) 280} 281 282 283qqmath.mcmc.list <- 284 function(x, data = NULL, 285 outer = FALSE, groups = !outer, 286 aspect = "xy", 287 default.scales = list(y = list(relation = "free")), 288 prepanel = prepanel.qqmathline, 289 start = 1, thin = 1, 290 main = attr(x, "title"), 291 ylab = "", 292 ..., 293 subset = thinned.indices(x[[1]], start = start, thin = thin)) 294{ 295 if (!is.R()) { 296 stop("This function is not yet available in S-PLUS") 297 } 298 if (groups && outer) warning("'groups=TRUE' ignored when 'outer=TRUE'") 299 datalist <- lapply(x, function(x) as.data.frame(x)[subset, , drop=FALSE]) 300 data <- do.call("rbind", datalist) 301 form <- 302 if (outer) 303 as.formula(paste("~", 304 paste(lapply(names(data), as.name), 305 collapse = "+"), 306 "| .run")) 307 else 308 as.formula(paste("~", 309 paste(lapply(names(data), as.name), 310 collapse = "+"))) 311 ##data[["index"]] <- seq(length = nrow(x[[1]]))[subset] 312 ##data[[".run"]] <- gl(length(datalist), nrow(datalist[[1]])) 313 .run <- gl(length(datalist), nrow(datalist[[1]])) 314 if (groups && !outer) 315 qqmath(form, data = data, 316 outer = TRUE, 317 groups = .run, 318 aspect = aspect, 319 prepanel = prepanel, 320 default.scales = default.scales, 321 main = main, 322 ylab = ylab, 323 ...) 324 else 325 qqmath(form, data = data, 326 outer = TRUE, 327 aspect = aspect, 328 default.scales = default.scales, 329 main = main, 330 ylab = ylab, 331 ...) 332} 333 334 335 336### methods for xyplot (mcmc and mcmc.list) 337 338 339 340xyplot.mcmc <- 341 function(x, data = NULL, 342 outer, layout = c(1, nvar(x)), 343 default.scales = list(y = list(relation = "free")), 344 type = 'l', 345 start = 1, thin = 1, 346 xlab = "Iteration number", 347 ylab = "", 348 main = attr(x, "title"), 349 ..., 350 subset = thinned.indices(x, start = start, thin = thin)) 351{ 352 if (!is.R()) { 353 stop("This function is not yet available in S-PLUS") 354 } 355 if (!missing(outer)) warning("specification of outer ignored") 356 data <- as.data.frame(x) 357 form <- eval(parse(text = paste(paste(lapply(names(data), as.name), 358 collapse = "+"), "~.index"))) 359 data[[".index"]] <- time(x) 360 xyplot(form, data = data[subset, ], 361 outer = TRUE, 362 layout = layout, 363 default.scales = default.scales, 364 type = type, 365 xlab = xlab, 366 ylab = ylab, 367 main = main, 368 ...) 369} 370 371 372 373xyplot.mcmc.list <- 374 function(x, data = NULL, 375 outer = FALSE, groups = !outer, 376 aspect = "xy", layout = c(1, nvar(x)), 377 default.scales = list(y = list(relation = "free")), 378 type = 'l', 379 start = 1, thin = 1, 380 xlab = "Iteration number", 381 ylab = "", 382 main = attr(x, "title"), 383 ..., 384 subset = thinned.indices(x[[1]], start = start, thin = thin)) 385{ 386 if (!is.R()) { 387 stop("This function is not yet available in S-PLUS") 388 } 389 if (groups && outer) warning("'groups=TRUE' ignored when 'outer=TRUE'") 390 datalist <- lapply(x, function(x) as.data.frame(x)[subset,,drop=FALSE]) 391 data <- do.call("rbind", datalist) 392 form <- 393 if (outer) 394 eval(parse(text = paste(paste(lapply(names(data), as.name), 395 collapse = "+"), "~.index | .run"))) 396 else 397 eval(parse(text = paste(paste(lapply(names(data), as.name), 398 collapse = "+"), "~.index"))) 399## form <- 400## if (outer) 401## as.formula(paste(paste(names(data), 402## collapse = "+"), 403## "~ index | .run")) 404## else 405## as.formula(paste(paste(names(data), 406## collapse = "+"), 407## "~ index")) 408 data[[".index"]] <- time(x) 409 .run <- gl(length(datalist), nrow(datalist[[1]])) 410 if (groups && !outer) 411 xyplot(form, data = data, 412 outer = TRUE, 413 layout = layout, 414 groups = .run, 415 default.scales = default.scales, 416 type = type, 417 main = main, 418 xlab = xlab, 419 ylab = ylab, 420 ...) 421 else 422 xyplot(form, data = data, 423 outer = TRUE, 424 layout = layout, 425 default.scales = default.scales, 426 type = type, 427 main = main, 428 xlab = xlab, 429 ylab = ylab, 430 ...) 431} 432 433 434 435 436 437 438 439 440### methods for acfplot (mcmc and mcmc.list) 441 442 443 444panel.acfplot <- 445 function(..., groups = NULL) 446{ 447 if (!is.R()) { 448 stop("This function is not yet available in S-PLUS") 449 } 450 reference.line <- trellis.par.get("reference.line") 451 panel.abline(h = 0, 452 col = reference.line$col, 453 lty = reference.line$lty, 454 lwd = reference.line$lwd, 455 alpha = reference.line$alpha) 456 if (is.null(groups)) panel.xyplot(...) 457 else panel.superpose(..., groups = groups) 458} 459 460 461 462 463acfplot <- function(x, data, ...) 464 UseMethod("acfplot") 465 466 467 468acfplot.mcmc <- 469 function(x, data = NULL, 470 outer, 471 prepanel = function(x, y, ...) list(ylim= c(-1, 1) * max(abs(y[-1]))), 472 panel = panel.acfplot, 473 type = "h", 474 aspect = "xy", 475 start = 1, thin = 1, 476 lag.max = NULL, 477 ylab = "Autocorrelation", 478 xlab = "Lag", 479 main = attr(x, "title"), 480 ..., 481 subset = thinned.indices(x, start = start, thin = thin)) 482{ 483 if (!is.R()) { 484 stop("This function is not yet available in S-PLUS") 485 } 486 if (!missing(outer)) warning("specification of outer ignored") 487 getAcf <- function(x, lag.max) 488 { 489 as.vector(acf(x, lag.max = lag.max, plot = FALSE)$acf) 490 } 491 data <- as.data.frame(apply(as.matrix(x)[subset, ,drop=FALSE], 2, getAcf, lag.max = lag.max)) 492 form <- 493 eval(parse(text = paste(paste(lapply(names(data), as.name), 494 collapse = "+"), "~.lag"))) 495 data[[".lag"]] <- seq(length = nrow(data)) 496 xyplot(form, data = data, 497 outer = TRUE, 498 prepanel = prepanel, 499 panel = panel, 500 type = type, 501 aspect = aspect, 502 xlab = xlab, 503 ylab = ylab, 504 main = main, 505 ...) 506} 507 508 509 510 511acfplot.mcmc.list <- 512 function(x, data = NULL, 513 outer = FALSE, groups = !outer, 514 prepanel = function(x, y, ..., groups = NULL, subscripts) { 515 if (is.null(groups)) list(ylim= c(-1, 1) * max(abs(y[-1]))) 516 else list(ylim = c(-1, 1) * max(sapply(split(y, groups[subscripts]), 517 function(x) 518 max(abs(x[-1]), na.rm = TRUE )))) 519 }, 520 panel = panel.acfplot, 521 type = if (groups) 'b' else 'h', 522 aspect = "xy", 523 start = 1, thin = 1, 524 lag.max = NULL, 525 ylab = "Autocorrelation", 526 xlab = "Lag", 527 main = attr(x, "title"), 528 ..., 529 subset = thinned.indices(x[[1]], start = start, thin = thin)) 530{ 531 if (!is.R()) { 532 stop("This function is not yet available in S-PLUS") 533 } 534 if (groups && outer) warning("'groups=TRUE' ignored when 'outer=TRUE'") 535 getAcf <- function(x, lag.max) 536 { 537 as.vector(acf(x, lag.max = lag.max, plot = FALSE)$acf) 538 } 539 if (groups || outer) 540 { 541 datalist <- 542 lapply(x, function(x) 543 as.data.frame(apply(as.matrix(x)[subset, ,drop=FALSE], 2, 544 getAcf, lag.max = lag.max))) 545 data <- do.call("rbind", datalist) 546 } 547 else 548 { 549 ## this is not quite valid, as we are combining multiple 550 ## series, but shouldn't be too bad (FIXME: should we warn?) 551 552 datalist <- lapply(x, function(x) as.matrix(x)[subset, ,drop=FALSE]) 553 data <- 554 as.data.frame(apply(do.call("rbind", datalist), 555 2, getAcf, lag.max = lag.max)) 556 } 557 form <- 558 if (outer) 559 as.formula(paste(paste(lapply(names(data), as.name), 560 collapse = "+"), 561 "~ .lag | .run")) 562 else 563 as.formula(paste(paste(lapply(names(data), as.name), 564 collapse = "+"), 565 "~ .lag")) 566 data[[".lag"]] <- seq(length = nrow(datalist[[1]])) ## repeated 567 .run <- gl(length(datalist), nrow(datalist[[1]])) 568 if (groups && !outer) 569 xyplot(form, data = data, 570 outer = TRUE, 571 groups = .run, 572 prepanel = prepanel, 573 panel = panel, 574 type = type, 575 aspect = aspect, 576 xlab = xlab, 577 ylab = ylab, 578 main = main, 579 ...) 580 else 581 xyplot(form, data = data, 582 outer = TRUE, 583 prepanel = prepanel, 584 panel = panel, 585 type = type, 586 aspect = aspect, 587 xlab = xlab, 588 ylab = ylab, 589 main = main, 590 ...) 591} 592 593 594