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