## ----------------------
## Banana shaped distribution
## log-pdf to sample from
p.log <- function(x) {
B <- 0.03 # controls 'bananacity'
-x[1]^2/200 - 1/2*(x[2]+B*x[1]^2-100*B)^2
}
## ----------------------
## generate samples
## 1) non-adaptive sampling
samp.1 <- MCMC(p.log, n=200, init=c(0, 1), scale=c(1, 0.1),
adapt=FALSE)
## 2) adaptive sampling
samp.2 <- MCMC(p.log, n=200, init=c(0, 1), scale=c(1, 0.1),
adapt=TRUE, acc.rate=0.234)
## ----------------------
## summarize results
str(samp.2)
summary(samp.2$samples)
## covariance of last jump distribution
samp.2$cov.jump
## ----------------------
## plot density and samples
x1 <- seq(-15, 15, length=80)
x2 <- seq(-15, 15, length=80)
d.banana <- matrix(apply(expand.grid(x1, x2), 1, p.log), nrow=80)
par(mfrow=c(1,2))
image(x1, x2, exp(d.banana), col=cm.colors(60), asp=1, main="no adaption")
contour(x1, x2, exp(d.banana), add=TRUE, col=gray(0.6))
lines(samp.1$samples, type='b', pch=3)
image(x1, x2, exp(d.banana), col=cm.colors(60), asp=1, main="with adaption")
contour(x1, x2, exp(d.banana), add=TRUE, col=gray(0.6))
lines(samp.2$samples, type='b', pch=3)
## ----------------------
## function returning extra information in a list
p.log.list <- function(x) {
B <- 0.03 # controls 'bananacity'
log.density <- -x[1]^2/200 - 1/2*(x[2]+B*x[1]^2-100*B)^2
result <- list(log.density=log.density)
## under some conditions one may want to return other infos
if(x[1]<0) {
result$message <- "Attention x[1] is negative!"
result$x <- x[1]
}
result
}
samp.list <- MCMC(p.log.list, n=200, init=c(0, 1), scale=c(1, 0.1),
adapt=TRUE, acc.rate=0.234)
## the additional values are stored under `extras.values`
head(samp.list$extras.values)
Run the code above in your browser using DataLab