foo <- read.table(url("http://www.stat.umn.edu/geyer/PhD/F03/logit.txt"),
    header = TRUE)

x <- foo
x$y <- NULL
x <- as.matrix(x)
x <- cbind(1, x)
dimnames(x) <- NULL
y <- foo$y

lupost <- function(beta) {
    eta <- x %*% beta
    p <- 1 / (1 + exp(- eta))   # note: works even when eta is Inf or -Inf
    logl <- sum(log(p[y == 1])) + sum(log(1 - p[y == 0]))
    return(logl + sum(dnorm(beta, 0, 2, log = TRUE)))
}

library(mcmc)
set.seed(42)
beta.init <- rep(0, ncol(x))
out <- metrop(lupost, beta.init, blen = 100, nbatch = 1e2, scale = 0.4)
out$accept

fred <- function(z) {
    sally <- 1 / (1 + exp(- z[1]))
    c(sally, sally^2)
}

out <- metrop(out, nbatch = 1e3, outfun = fred)
out$accept

# plot(out$batch[ , 1])
# plot(out$batch[ , 2])
# acf(out$batch[ , 1])
# acf(out$batch[ , 2])
# barely signif lag-1, double batch length?

mu.hat <- mean(out$batch[ , 1])
mu.mcse <- sd(out$batch[ , 1]) / sqrt(out$nbatch)
print(mu.hat)
print(mu.mcse)

var.hat <- mean(out$batch[ , 2] - mu.hat^2)
##### bogus
var.mcse <- sd(out$batch[ , 2] - mu.hat^2) / sqrt(out$nbatch)
##### right
delta1 <- out$batch[ , 1] - mean(out$batch[ , 1])
delta2 <- out$batch[ , 2] - mean(out$batch[ , 2])
var.mcse.ok <- sqrt(mean((delta2 - 2 * mu.hat * delta1)^2) / out$nbatch)
print(var.hat)
print(var.mcse)
print(var.mcse.ok)

sigma.hat <- sqrt(var.hat)
sigma.mcse <- var.mcse / (2 * sigma.hat)
sigma.mcse.ok <- var.mcse.ok / (2 * sigma.hat)
print(sigma.hat)
print(sigma.mcse)
print(sigma.mcse.ok)

##### looks like merely doubling the batch length should do it

out <- metrop(out, blen = 250, outfun = fred)
out$accept

mu.hat <- mean(out$batch[ , 1])
mu.mcse <- sd(out$batch[ , 1]) / sqrt(out$nbatch)
print(mu.hat)
print(mu.mcse)

var.hat <- mean(out$batch[ , 2] - mu.hat^2)
##### bogus
var.mcse <- sd(out$batch[ , 2] - mu.hat^2) / sqrt(out$nbatch)
##### right
delta1 <- out$batch[ , 1] - mean(out$batch[ , 1])
delta2 <- out$batch[ , 2] - mean(out$batch[ , 2])
var.mcse.ok <- sqrt(mean((delta2 - 2 * mu.hat * delta1)^2) / out$nbatch)
print(var.hat)
print(var.mcse)
print(var.mcse.ok)

sigma.hat <- sqrt(var.hat)
sigma.mcse <- var.mcse / (2 * sigma.hat)
sigma.mcse.ok <- var.mcse.ok / (2 * sigma.hat)
print(sigma.hat)
print(sigma.mcse)
print(sigma.mcse.ok)

print(out$nbatch)
print(out$blen)
print(out$nspac)
print(out$nburn)
print(out$nbatch * out$blen)