MCMC
Janpu Hou
September 6, 2017
MCMC
n Simulates continuous distributions of random vectors using Markov chain Monte Carlo (MCMC). Users specify the distribution by an R function that evaluates the log unnormalized density. Algorithms are random walk Metropolis algorithm (function metrop), simulated tempering (function temper), and morphometric random walk Metropolis
library(mcmc)
data(foo)
out <- glm(y ~ x1 + x2 + x3, family = binomial, data = foo)
summary(out)##
## Call:
## glm(formula = y ~ x1 + x2 + x3, family = binomial, data = foo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0371 -0.6337 0.2394 0.6685 1.9599
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5772 0.2766 2.087 0.036930 *
## x1 0.3362 0.4256 0.790 0.429672
## x2 0.8475 0.4701 1.803 0.071394 .
## x3 1.5143 0.4426 3.422 0.000622 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 134.602 on 99 degrees of freedom
## Residual deviance: 86.439 on 96 degrees of freedom
## AIC: 94.439
##
## Number of Fisher Scoring iterations: 5Initial Sequence Estimators
Variance of sample mean of functional of reversible Markov chain using methods of Geyer (1992).
n <- 2e4
rho <- 0.99
x <- arima.sim(model = list(ar = rho), n = n)
out <- initseq(x)
## Not run:
plot(seq(along = out$Gamma.pos) - 1, out$Gamma.pos,
xlab = "k", ylab = expression(Gamma[k]), type = "l")
lines(seq(along = out$Gamma.dec) - 1, out$Gamma.dec, col = "red")
lines(seq(along = out$Gamma.con) - 1, out$Gamma.con, col = "blue")
## End(Not run)
# asymptotic 95% confidence interval for mean of x
mean(x) + c(-1, 1) * qnorm(0.975) * sqrt(out$var.con / length(x))## [1] -1.059428 1.364819# estimated asymptotic variance
out$var.con## [1] 7649.409# theoretical asymptotic variance
(1 + rho) / (1 - rho) * 1 / (1 - rho^2)## [1] 10000# illustrating use with batch means
bm <- apply(matrix(x, nrow = 5), 2, mean)
initseq(bm)$var.con * 5## [1] 7663.986logit Simulated logistic regression data.
data(logit)
out <- glm(y ~ x1 + x2 + x3 + x4, family = binomial, data = logit)
summary(out)##
## Call:
## glm(formula = y ~ x1 + x2 + x3 + x4, family = binomial, data = logit)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7461 -0.6907 0.1540 0.7041 2.1943
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6328 0.3007 2.104 0.03536 *
## x1 0.7390 0.3616 2.043 0.04100 *
## x2 1.1137 0.3627 3.071 0.00213 **
## x3 0.4781 0.3538 1.351 0.17663
## x4 0.6944 0.3989 1.741 0.08172 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 137.628 on 99 degrees of freedom
## Residual deviance: 87.668 on 95 degrees of freedom
## AIC: 97.668
##
## Number of Fisher Scoring iterations: 6metrop Metropolis Algorithm
Markov chain Monte Carlo for continuous random vector using a Metropolis algorithm.
h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
out <- metrop(h, rep(0, 5), 1000)
out$accept## [1] 0# acceptance rate too low
out <- metrop(out, scale = 0.1)
out$accept## [1] 0.174t.test(out$accept.batch)$conf.int## [1] 0.1504627 0.1975373
## attr(,"conf.level")
## [1] 0.95# acceptance rate o. k. (about 25 percent)
plot(out$batch[ , 1])
# but run length too short (few excursions from end to end of range)
out <- metrop(out, nbatch = 1e4)
out$accept## [1] 0.2246plot(out$batch[ , 1])
hist(out$batch[ , 1])
acf(out$batch[ , 1], lag.max = 250)
# looks like batch length of 250 is perhaps OK
out <- metrop(out, blen = 250, nbatch = 100)
apply(out$batch, 2, mean) # Monte Carlo estimates of means## [1] 0.1702904 0.1691711 0.1637148 0.1687371 0.1625102apply(out$batch, 2, sd) / sqrt(out$nbatch) # Monte Carlo standard errors## [1] 0.006773083 0.007693173 0.006337058 0.006863295 0.005705977t.test(out$accept.batch)$conf.int## [1] 0.2166265 0.2318535
## attr(,"conf.level")
## [1] 0.95acf(out$batch[ , 1]) # appears that blen is long enough
morph Variable Transformation
# use an exponential transformation, centered at 100.
b1 <- morph(b=1, center=100)
# original log unnormalized density is from a t distribution with 3
# degrees of freedom, centered at 100.
lud.transformed <- b1$lud(function(x) dt(x - 100, df=3, log=TRUE))
d.transformed <- Vectorize(function(x) exp(lud.transformed(x)))
## Not run:
curve(d.transformed, from=-3, to=3, ylab="Induced Density")
## End(Not run)morph.metrop Morphometric Metropolis Algorithm
Markov chain Monte Carlo for continuous random vector using a Metropolis algorithm for an induced density.
out <- morph.metrop(function(x) dt(x, df=3, log=TRUE), 0, blen=100,
nbatch=100, morph=morph(b=1))
# change the transformation.
out <- morph.metrop(out, morph=morph(b=2))
out$accept## [1] 0.6291# accept rate is high, increase the scale.
out <- morph.metrop(out, scale=4)
# close to 0.20 is about right.
out$accept## [1] 0.2238olbm Overlapping Batch Means
Variance of sample mean of time series calculated using overlapping batch means
h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
out <- metrop(h, rep(0, 5), 1000)
out <- metrop(out, scale = 0.1)
out <- metrop(out, nbatch = 1e4)
foo <- olbm(out$batch, 150)
# monte carlo estimates (true means are same by symmetry)
apply(out$batch, 2, mean)## [1] 0.1660160 0.1575581 0.1725433 0.1603277 0.1743106# monte carlo standard errors (true s. d. are same by symmetry)
sqrt(diag(foo))## [1] 0.009815753 0.008866610 0.011476625 0.009300513 0.012614992# check that batch length is reasonable
acf(out$batch, lag.max = 200)
temper Simulated Tempering and Umbrella Sampling
Markov chain Monte Carlo (MCMC) for continuous random vectors using parallel or serial tempering, the latter also called umbrella sampling and simulated tempering. The chain simulates k different distributions on the same state space. In parallel tempering, all the distributions are simulated in each iteration. In serial tempering, only one of the distributions is simulated (a random one). In parallel tempering, the state is a k × p matrix, each row of which is the state for one of the distributions. In serial tempering the state of the Markov chain is a pair (i, x), where i is an integer between 1 and k and x is a vector of length p. This pair is represented as a single real vector c(i, x). The variable i says which distribution x is a simulation for.
d <- 9
witch.which <- c(0.1, 0.3, 0.5, 0.7, 1.0)
ncomp <- length(witch.which)
neighbors <- matrix(FALSE, ncomp, ncomp)
neighbors[row(neighbors) == col(neighbors) + 1] <- TRUE
neighbors[row(neighbors) == col(neighbors) - 1] <- TRUE
ludfun <- function(state, log.pseudo.prior = rep(0, ncomp)) {
stopifnot(is.numeric(state))
stopifnot(length(state) == d + 1)
icomp <- state[1]
stopifnot(icomp == as.integer(icomp))
stopifnot(1 <= icomp && icomp <= ncomp)
stopifnot(is.numeric(log.pseudo.prior))
stopifnot(length(log.pseudo.prior) == ncomp)
theta <- state[-1]
if (any(theta > 1.0)) return(-Inf)
bnd <- witch.which[icomp]
lpp <- log.pseudo.prior[icomp]
if (any(theta > bnd)) return(lpp)
return(- d * log(bnd) + lpp)
}
# parallel tempering
thetas <- matrix(0.5, ncomp, d)
out <- temper(ludfun, initial = thetas, neighbors = neighbors, nbatch = 20,
blen = 10, nspac = 5, scale = 0.56789, parallel = TRUE, debug = TRUE)
# serial tempering
theta.initial <- c(1, rep(0.5, d))
# log pseudo prior found by trial and error
qux <- c(0, 9.179, 13.73, 16.71, 20.56)
out <- temper(ludfun, initial = theta.initial, neighbors = neighbors,
nbatch = 50, blen = 30, nspac = 2, scale = 0.56789,
parallel = FALSE, debug = FALSE, log.pseudo.prior = qux)
plot(out$batch[ , 1])
hist(out$batch[ , 1])
acf(out$batch[ , 1], lag.max = 250)