R/heidel.R In coda: Output Analysis and Diagnostics for MCMC

```"heidel.diag" <- function (x, eps = 0.1, pvalue=0.05)
{
if (is.mcmc.list(x))
return(lapply(x, heidel.diag, eps))
x <- as.mcmc(as.matrix(x))
HW.mat0 <- matrix(0, ncol = 6, nrow = nvar(x))
dimnames(HW.mat0) <- list(varnames(x),
c("stest", "start", "pvalue", "htest",
"mean", "halfwidth"))
HW.mat <- HW.mat0
for (j in 1:nvar(x)) {
start.vec <- seq(from=start(x), to = end(x)/2, by=niter(x)/10)
Y <- x[, j, drop = TRUE]
n1 <- length(Y)
## Schruben's test for convergence, applied sequentially
##
S0 <- spectrum0.ar(window(Y, start=end(Y)/2))\$spec
converged <- FALSE
for (i in seq(along = start.vec)) {
Y <- window(Y, start = start.vec[i])
n <- niter(Y)
ybar <- mean(Y)
B <- cumsum(Y) - ybar * (1:n)
Bsq <- (B * B)/(n * S0)
I <- sum(Bsq)/n
if(converged <- !is.na(I) && pcramer(I) < 1 - pvalue)
break
}
## Recalculate S0 using section of chain that passed convergence test
S0ci <- spectrum0.ar(Y)\$spec
halfwidth <- 1.96 * sqrt(S0ci/n)
passed.hw <- !is.na(halfwidth) & (abs(halfwidth/ybar) <= eps)
if (!converged || is.na(I) || is.na(halfwidth)) {
nstart <- NA
passed.hw <- NA
halfwidth <- NA
ybar <- NA
}
else {
nstart <- start(Y)
}
HW.mat[j, ] <- c(converged, nstart, 1 - pcramer(I),
passed.hw, ybar, halfwidth)
}
class(HW.mat) <- "heidel.diag"
return(HW.mat)
}

"print.heidel.diag" <-
function (x, digits = 3, ...)
{
HW.title <- matrix(c("Stationarity", "test", "start", "iteration",
"p-value", "",
"Halfwidth", "test", "Mean", "", "Halfwidth", ""),
nrow = 2)
y <- matrix("", nrow = nrow(x), ncol = 6)
for (j in 1:ncol(y)) {
y[, j] <- format(x[, j], digits = digits)
}
y[, c(1, 4)] <- ifelse(x[, c(1, 4)], "passed", "failed")
y <- rbind(HW.title, y)
vnames <- if (is.null(rownames(x)))
paste("[,", 1:nrow(x), "]", sep = "")
else rownames(x)
dimnames(y) <- list(c("", "", vnames), rep("", 6))
print.default(y[, 1:3], quote = FALSE, ...)
print.default(y[, 4:6], quote = FALSE, ...)
invisible(x)
}

"spectrum0.ar" <- function(x)
{
x <- as.matrix(x)
v0 <- order <- numeric(ncol(x))
names(v0) <- names(order) <- colnames(x)
z <- 1:nrow(x)
for (i in 1:ncol(x))
{
lm.out <- lm(x[,i] ~ z)
if (identical(all.equal(sd(residuals(lm.out)), 0), TRUE)) {
v0[i] <- 0
order[i] <- 0
}
else {
ar.out <- ar(x[,i], aic=TRUE)
v0[i] <- ar.out\$var.pred/(1 - sum(ar.out\$ar))^2
order[i] <- ar.out\$order
}
}
return(list(spec=v0, order=order))
}

effectiveSize <- function(x)
{
if (is.mcmc.list(x))
{
##RGA changed to sum across all chains
ess <- do.call("rbind",lapply(x,effectiveSize))
ans <- apply(ess,2,sum)
}
else
{
x <- as.mcmc(x)
x <- as.matrix(x)
spec <- spectrum0.ar(x)\$spec
ans <- ifelse(spec==0, 0, nrow(x) * apply(x, 2, var)/spec)
}
return(ans)
}

"spectrum0" <- function(x, max.freq=0.5, order=1, max.length=200)
{
x <- as.matrix(x)
if (!is.null(max.length) && nrow(x) > max.length) {
batch.size <- ceiling(nrow(x)/max.length)
if (is.R()) {
x <- aggregate(ts(x, frequency=batch.size), nfreq = 1, FUN=mean)
}
else {
x <- aggregate(ts(x, frequency=batch.size), nf = 1, fun=mean)
}
}
else {
batch.size <- 1
}

out <- do.spectrum0(x, max.freq=max.freq, order=order)
out\$spec <- out\$spec * batch.size
return(out)
}

"do.spectrum0" <- function(x, max.freq=0.5, order=1)
{
## Estimate spectral density of time series x at frequency 0.
## spectrum0(x)/length(x) estimates the variance of mean(x)
##
## NB We do NOT use the same definition of spectral density
## as in spec.pgram.
##
fmla <- switch(order+1,
spec ~ one,
spec ~ f1,
spec ~ f1 + f2)
if(is.null(fmla))
stop("invalid order")

N <- nrow(x)
Nfreq <- floor(N/2)
freq <- seq(from = 1/N, by = 1/N, length = Nfreq)
f1 <- sqrt(3) * (4 * freq - 1)
f2 <- sqrt(5) * (24 * freq^2 - 12 * freq + 1)
v0 <- numeric(ncol(x))
for(i in 1:ncol(x)) {
y <- x[,i]
if (var(y) == 0) {
v0[i] <- 0
}
else {
yfft <- fft(y)
spec <- Re(yfft * Conj(yfft))/ N
spec.data <- data.frame(one = rep(1, Nfreq), f1=f1, f2=f2,
spec = spec[1 + (1:Nfreq)],
inset = I(freq<=max.freq))

v0[i] <- predict(glm.out, type="response",
newdata=data.frame(spec=0,one=1,f1=-sqrt(3),f2=sqrt(5)))
}
}
return(list(spec=v0))
}

"pcramer" <- function (q, eps=1.0e-5)
{
## Distribution function of the Cramer-von Mises statistic
##
log.eps <- log(eps)
y <- matrix(0, nrow=4, ncol=length(q))
for(k in 0:3) {
z <- gamma(k + 0.5) * sqrt(4*k + 1)/(gamma(k+1) * pi^(3/2) * sqrt(q))
u <- (4*k + 1)^2/(16*q)
y[k+1,] <- ifelse(u > -log.eps, 0, z * exp(-u) * besselK(x = u, nu=1/4))
}
return(apply(y,2,sum))
}
```

Try the coda package in your browser

Any scripts or data that you put into this service are public.

coda documentation built on July 5, 2019, 5:03 p.m.