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