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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
"geweke.diag" <-
function (x, frac1 = 0.1, frac2 = 0.5)
##
{
if (is.mcmc.list(x))
return(lapply(x, geweke.diag, frac1, frac2))
x <- as.mcmc(x)
xstart <- c(start(x), end(x) - frac2 * (end(x) - start(x)))
xend <- c(start(x) + frac1 * (end(x) - start(x)), end(x))
y.variance <- y.mean <- vector("list", 2)
for (i in 1:2) {
y <- window(x, start = xstart[i], end = xend[i])
y.mean[[i]] <- apply(as.matrix(y), 2, mean)
y.variance[[i]] <- spectrum0(y)$spec/niter(y)
}
z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])
out <- list(z = z, frac = c(frac1, frac2))
class(out) <- "geweke.diag"
return(out)
}
"geweke.plot" <-
function (x, frac1 = 0.1, frac2 = 0.5, nbins = 20,
pvalue = 0.05, auto.layout = TRUE, ask, ...)
{
if (missing(ask)) {
ask <- if (is.R()) {
dev.interactive()
}
else {
interactive()
}
}
x <- as.mcmc.list(x)
oldpar <- NULL
on.exit(par(oldpar))
if (auto.layout)
oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x),
Nparms = nvar(x)))
ystart <- seq(from = start(x), to = (start(x) + end(x))/2, length = nbins)
if (is.R())
gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)),
dimnames = c(ystart, varnames(x), chanames(x)))
else
gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)),
dimnames = list(ystart, varnames(x), chanames(x)))
for (n in 1:length(ystart)) {
geweke.out <- geweke.diag(window(x, start = ystart[n]),
frac1 = frac1, frac2 = frac2)
for (k in 1:nchain(x)) gcd[n, , k] <- geweke.out[[k]]$z
}
climit <- qnorm(1 - pvalue/2)
for (k in 1:nchain(x)) for (j in 1:nvar(x)) {
ylimit <- max(c(climit, abs(gcd[, j, k])))
plot(ystart, gcd[, j, k], type = "p", xlab = "First iteration in segment",
ylab = "Z-score", pch = 4, ylim = c(-ylimit, ylimit),
...)
abline(h = c(climit, -climit), lty = 2)
if (nchain(x) > 1) {
title(main = paste(varnames(x, allow.null = FALSE)[j],
" (", chanames(x, allow.null = FALSE)[k], ")",
sep = ""))
}
else {
title(main = paste(varnames(x, allow.null = FALSE)[j],
sep = ""))
}
if (k==1 && j==1)
oldpar <- c(oldpar, par(ask = ask))
}
invisible(list(start.iter = ystart, z = gcd))
}
"print.geweke.diag" <- function (x, digits = min(4, .Options$digits), ...)
## Print method for output from geweke.diag
{
cat("\nFraction in 1st window =", x$frac[1])
cat("\nFraction in 2nd window =", x$frac[2], "\n\n")
print.default(x$z, digits = digits, ...)
cat("\n")
invisible(x)
}
|
